🔢 mathematics

GPU Performance of an Entropy-Stable Discontinuous Galerkin Euler Solver with Non-Conservative Terms

This paper presents a highly optimized GPU implementation of an entropy-stable discontinuous Galerkin Euler solver with non-conservative buoyancy terms, demonstrating a 10x speedup and over 13x improvement in energy efficiency compared to a CPU code while achieving nearly 70% of the A100's double-precision peak performance.

Original authors: Henry Waterhouse, Maciej Waruszewski, Lucas C. Wilcox, Francis X. Giraldo

Published 2026-05-19

📖 4 min read🧠 Deep dive

CC BY 4.0

Original authors: Henry Waterhouse, Maciej Waruszewski, Lucas C. Wilcox, Francis X. Giraldo

Original paper licensed under CC BY 4.0 (http://creativecommons.org/licenses/by/4.0/). ✨ This is an AI-generated explanation of the paper below. It is not written or endorsed by the authors. For technical accuracy, refer to the original paper. Read full disclaimer

Imagine you are trying to predict the weather. To do this, scientists use complex computer models that simulate how air moves, heats up, and cools down. These simulations are like massive, digital wind tunnels. However, running these simulations accurately is incredibly difficult and slow, often requiring supercomputers that consume huge amounts of electricity.

This paper is about a team of researchers who built a new, super-fast version of one of these weather simulation tools. They took a specific mathematical method called "Entropy-Stable Discontinuous Galerkin" (a fancy name for a very precise way of calculating air movement) and rewrote it to run on GPUs (Graphics Processing Units).

Think of a CPU (the brain of a standard computer) as a single, very smart chef who can cook a gourmet meal perfectly but can only do one thing at a time. A GPU, on the other hand, is like a kitchen with thousands of line cooks. Each cook is less "smart" individually, but together, they can chop, fry, and plate thousands of ingredients simultaneously.

Here is what the researchers achieved, broken down into simple concepts:

1. The Problem: The "Two-Point" Bottleneck

The mathematical method they used is great at keeping the simulation stable (it doesn't crash or produce nonsense), but it's usually very slow. It's like trying to calculate the weather by asking every single air molecule to talk to every other molecule next to it. This involves a lot of heavy math, specifically "two-point fluxes," which are like complex handshake calculations between neighbors.

On a standard computer, this process is so heavy that it slows everything down. The researchers wanted to see if they could make this "handshake" process fast enough to run on the thousands of cooks in a GPU kitchen.

2. The Solution: Optimizing the GPU Kitchen

The team didn't just move the code to the GPU; they completely reorganized the kitchen to make it efficient. They made several clever changes:

Pre-cooking Ingredients: Some math operations, like logarithms and divisions, are very slow on GPUs. The researchers realized they could calculate these once and save them, rather than recalculating them every single time. It's like pre-chopping all the onions before the cooking starts instead of chopping one onion every time you need a pinch.
Symmetry Tricks: Many of the calculations are mirror images of each other. Instead of doing the math for both sides, they figured out how to do it for one side and just flip the result. This cut their work in half.
Balancing the Workload: In a kitchen with 1,000 cooks, you don't want 500 cooks doing nothing while the other 500 are overwhelmed. They created a system to ensure every "cook" (GPU thread) had exactly the same amount of work to do, preventing bottlenecks.

3. The Results: Speed and Efficiency

When they tested their new system on a powerful NVIDIA A100 GPU (a top-tier scientific chip), the results were impressive:

Raw Speed: The new GPU code ran 10 times faster than their highly optimized code running on a standard CPU.
Energy Efficiency: Because it was so much faster, it used 13 times less energy to do the same job. This is like getting 13 miles per gallon instead of 1 mile per gallon.
Peak Performance: Their system was able to use nearly 70% of the GPU's maximum theoretical speed, which is a very high score for such complex math.

4. Real-World Tests

They didn't just run numbers; they simulated two specific weather scenarios:

A Rising Thermal Bubble: Imagine a hot air balloon rising through cold air. They simulated this in 3D.
Baroclinic Instability: This is a complex weather pattern that leads to storms and fronts, like the ones that cause winter storms in the mid-latitudes.

They found that the GPU version produced results that were just as accurate as the CPU version. Interestingly, they also tested the simulation using "32-bit" math (a slightly less precise but faster way of calculating numbers). They found that while the results were almost identical to the high-precision version, the simulation ran twice as fast. This suggests that for many weather predictions, we might be able to use the faster, slightly less precise math without losing much accuracy.

Summary

In short, this paper shows that by treating a computer chip like a massive team of workers rather than a single genius, and by organizing their tasks with extreme efficiency, scientists can run weather simulations 10 times faster and use 13 times less electricity. This brings us closer to having weather models that are both highly accurate and fast enough to run on more accessible hardware.

Technical Summary: GPU Performance of an Entropy-Stable Discontinuous Galerkin Euler Solver with Non-Conservative Terms

Problem Statement

The construction of accurate, stable, and efficient dynamical cores for atmospheric modeling remains a critical challenge. While Discontinuous Galerkin (DG) methods offer high-order accuracy and excellent parallelization potential, their application to geophysical fluid dynamics faces two primary hurdles: stability at high orders and computational efficiency. Entropy-stable DG (ESDG) methods resolve stability issues by constructing numerical fluxes that conserve or dissipate entropy, but they are computationally expensive due to the requirement of two-point flux evaluations. Furthermore, incorporating buoyancy into the compressible Euler equations transforms the system from a set of conservation laws into a system of balance laws with non-conservative terms, complicating the flux formulation.

Despite extensive literature on GPU implementations for Euler equations and shallow water equations, there is a significant gap in published work regarding the GPU optimization of ESDG methods, particularly for the Euler equations with buoyancy. Without addressing the high computational cost, ESDG methods risk being impractical for large-scale atmospheric modeling. This paper addresses the need to construct efficient ESDG methods on modern Graphics Processing Unit (GPU) architectures.

Methodology

The authors implement an entropy-stable discontinuous Galerkin solver for the compressible Euler equations with buoyancy on NVIDIA A100 GPUs. The numerical approach follows the ESDG framework detailed in prior work [55], adapted for GPU execution.

Governing Equations and Discretization

The solver handles the time-dependent compressible Euler equations including geopotential (buoyancy) and Coriolis terms. The domain is discretized using non-overlapping hexahedral elements mapped to reference elements containing $(N+1)^3$ Legendre-Gauss-Lobatto (LGL) nodes. The semi-discrete form utilizes a matrix-free approach exploiting tensor-product structures to minimize memory loads. Entropy stability is enforced via:

Volume Terms: Using an entropy-conserving two-point flux.
Surface Terms: Using an entropy-conserving flux augmented with matrix dissipation (similar to a Roe Riemann solver).
Time Integration: A fourth-order, low-storage Runge-Kutta scheme.

GPU Optimization Strategies

The core of the paper focuses on a series of optimizations to mitigate the computational overhead of two-point fluxes and non-conservative terms:

Thread Mapping: Each element is assigned to a thread block, with threads representing degrees of freedom (DOFs).
Operation Reduction:
- Precomputation: Primitive variables (e.g., velocity from momentum) and logarithmic terms are precomputed once per thread rather than recalculated for every flux evaluation.
- Logarithmic Means: The authors adopt an efficient formulation [42] using Horner's rule to reduce divisions and rewrite logarithmic means to compute $\log(a)$ and $\log(b)$ individually, allowing reuse across evaluations.
Symmetry Exploitation:
- For standard Euler equations, symmetry in flux pairs allows a factor-of-two reduction in computations.
- Non-Symmetric Gravity Handling: The buoyancy term introduces non-symmetric components. The authors decompose the momentum flux into a symmetric part and a gravity contribution. The gravity term is computed once per pair and scaled for the partner node, preserving nearly the full factor-of-two speedup.
Load Balancing: To prevent load imbalance when threads compute flux pairs ( $j > i$ ), a specific indexing strategy is employed. For even polynomial orders, threads perform equal work; for odd orders, a weighting approach ensures uniform loop iterations, avoiding synchronization issues.
Memory and Kernel Fusion: State variables and precomputed data are stored in shared memory for fast access. The three spatial directions are fused into a single kernel launch to eliminate redundant shared memory loads. Contravariant fluxes are computed directly to avoid forming and contracting intermediate matrices.

Key Results

Performance on Single GPU

Compute Efficiency: On NVIDIA A100 hardware, the optimized volume kernel achieves nearly 70% of the 64-bit floating-point peak performance (roofline ceiling) for polynomial orders $N \ge 3$ . The kernel is compute-bound, not memory-bound, making it well-suited for GPU acceleration.
Optimization Impact: The cumulative optimizations (precomputation, efficient logarithmic means, and symmetry exploitation) yield a total speedup of 7.3× over the baseline GPU kernel.
Kernel Bottlenecks: Before optimization, the volume kernel dominated runtime. After optimization, the surface kernel becomes the dominant cost, indicating future optimization efforts should target memory-bound surface integrals.
Scaling: The solver demonstrates strong saturation capabilities. For high-order polynomials ( $N=6$ ), the GPU is fully saturated with as few as 4,000 elements.

Multi-GPU Scaling

Strong Scaling: Using the rising thermal bubble test case ( $N=4, L=6$ ), the solver scales up to 32 GPUs. Parallel efficiency drops to 60% at 32 GPUs due to communication overhead, suggesting optimal performance is achieved when problem size matches processor count.
Weak Scaling: When problem size and processor count increase proportionally, the solver maintains approximately constant time per step, demonstrating ideal weak scaling behavior.

GPU vs. CPU Comparison

The solver was compared against a highly optimized CPU code (NUMA model) running on AMD EPYC 7763 Milan nodes.

Speed: The GPU implementation is approximately 10× faster than the CPU code for the same problem size.
Energy Efficiency: The GPU implementation is >13× more energy efficient than the CPU code.
Precision Trade-offs:
- Rising Bubble: Switching from 64-bit to 32-bit precision yields a 2.16× speedup with minimal loss in large-scale symmetry.
- Baroclinic Instability: A smaller problem size resulted in a 1.5× speedup for 32-bit precision due to higher communication overhead relative to computation. However, increasing the problem size restored the 2.16× speedup. Visual results in 32-bit precision captured the same dynamics as 64-bit, with noise attributed to the magnitude of total energy variables rather than the method's stability.

Significance and Claims

The paper claims to present the first work on ESDG for the Euler equations focusing specifically on GPU optimization. Its primary significance lies in demonstrating that the computational cost of entropy-stable methods, often a barrier to their adoption, can be effectively mitigated on modern GPU hardware.

Key claims include:

Feasibility: ESDG methods for buoyancy-driven flows are viable on GPUs, achieving high arithmetic intensity and near-peak performance.
Optimization Roadmap: The paper provides a detailed roadmap of specific strategies (symmetry handling for non-conservative terms, load balancing, and operation reduction) that significantly improve efficiency.
Performance Gains: The implementation achieves a 10× speedup and >13× energy efficiency improvement over state-of-the-art CPU codes, making high-order, entropy-stable atmospheric modeling practical.
Precision Flexibility: The solver supports 32-bit precision with substantial speedups (up to 2×) while maintaining solution quality for geophysical problems, offering a practical trade-off for large-scale simulations.

The authors conclude that while the current work focuses on explicit time integration, the demonstrated spatial discretization performance lays the groundwork for future practical atmospheric models that may incorporate implicit time-stepping schemes.