A Massively Parallel Hybridizable Discontinuous Galerkin Solver for Direct Numerical Simulation of Compressible Flows on GPUs

Andrew Welter; Thea Collin; Ngoc Cuong Nguyen; Jaime Peraire

Source-linked AI summary

A Massively Parallel Hybridizable Discontinuous Galerkin Solver for Direct Numerical Simulation of Compressible Flows on GPUs

Andrew Welter, Thea Collin, Ngoc Cuong Nguyen, Jaime Peraire

arXiv:2608.27863v1Aug 28, 2026math.NAphysics.flu-dyn

TL;DR

DNS of compressible transitional and turbulent flows requires methods that jointly provide high-order accuracy, robustness, and computational efficiency across many scales and Mach regimes. The paper develops an integrated massively parallel HDG solver combining shock capturing, DIRK integration, Newton–GMRES, additive Schwarz preconditioning, reduced-basis initialization, and GPU execution. Across subsonic, supersonic, and hypersonic benchmarks, it reports agreement with experiments and published DNS data while resolving transition, compressibility, shock-associated structures, and three-dimensional turbulence.

Problem
Existing methods provide many relevant components separately, but an integrated massively parallel HDG methodology for repeated implicit DNS solves across subsonic, supersonic, and hypersonic flows remains insufficiently consolidated.
Method
The paper combines high-order HDG, localized shock capturing, DIRK time integration, Newton–GMRES, additive Schwarz preconditioning, reduced-basis initialization, and distributed GPU implementation.
Results
The solver agrees with experiments and published DNS data across three benchmark problems spanning subsonic, supersonic, and hypersonic regimes.
Takeaways & Limitations
The demonstrated framework resolves laminar–turbulent transition, strong compressibility effects, shock-associated structures, and fully three-dimensional turbulent dynamics.
Takeaways & Limitations
The hypersonic simulations use calorically perfect-gas equations; predictive high-enthalpy applications require real-gas thermochemistry and nonequilibrium extensions.

Abstract

from arXiv · show

Direct numerical simulation (DNS) of compressible transitional and turbulent flows requires numerical methods that combine high-order accuracy, robustness, and computational efficiency to resolve a broad range of spatial and temporal scales. This paper presents a massively parallel hybridizable discontinuous Galerkin (HDG) solver for DNS of the compressible Navier-Stokes equations on GPU-accelerated high-performance computing systems. The proposed solver combines high-order HDG discretization with robust shock capturing, diagonally implicit Runge-Kutta (DIRK) time integration, and an efficient Newton-GMRES solution strategy accelerated by additive Schwarz preconditioning and reduced-basis approximation. A distributed implementation of these methods based on GPU-aware MPI, Kokkos, and CUDA/HIP libraries enables scalable execution on heterogeneous computing platforms. The solver is demonstrated on three canonical benchmark problems covering a wide range of Mach-number flow regimes: subsonic transitional flow over the Eppler 387 airfoil, the supersonic Taylor-Green vortex, and hypersonic boundary-layer transition. Numerical results are compared with available experimental measurements and published DNS data, showing good agreement across distinct flow regimes. The results demonstrate the ability of the proposed solver to resolve laminar-turbulent transition, strong compressibility effects, shock-associated flow structures, and fully three-dimensional turbulent dynamics.

1. Introduction

The paper addresses the difficulty of performing accurate, robust, and scalable DNS of compressible transitional and turbulent flows across complex, shock-containing regimes. It integrates high-order HDG, localized shock treatment, implicit solvers, and GPU parallelism into one methodology.

DNS must resolve coupled spatial and temporal scales, instability growth, turbulence cascades, acoustic waves, and thermodynamic effects without excessive numerical dissipation.
High-order DG offers element-local approximation, compact coupling, conservation, low dispersion and dissipation, and suitability for complex meshes and parallel implementation.
Wall-resolved compressible DNS requires curved high-order elements and anisotropic meshes that concentrate degrees of freedom near walls, shocks, and rapidly varying solution regions.
Shock treatment must preserve high-order accuracy in smooth regions while adding dissipation locally to control oscillations near shocklets, compression waves, and strong shock structures.
Implicit DIRK integration avoids restrictive explicit stability limits but creates repeated nonlinear and linear systems whose cost can dominate fine-mesh DNS.
The proposed framework combines HDG, localized shock capturing, curved anisotropic meshes, DIRK, Newton–GMRES, additive Schwarz, reduced-basis initialization, and distributed GPU execution.

2. Numerical Methods

The numerical method discretizes the compressible Navier–Stokes equations with high-order HDG, localized artificial viscosity, DIRK time integration, and condensed trace-system solves. Element-local unknowns are eliminated to obtain a globally coupled trace problem that can be solved efficiently on distributed GPU platforms.

2.1. Governing Equations: The formulation solves the compressible Navier–Stokes equations in conservative form on a three-dimensional spatial domain without volumetric source terms.
2.1. Governing Equations: The conservative variables contain density, Cartesian velocity components, and total specific energy, with pressure, enthalpy, and transport relations defined for a calorically perfect gas.
2.1. Governing Equations: Artificial viscosity adds localized Laplacian diffusion to conservative variables and uses a compression sensor attenuated in predominantly vortical flow by a Ducros-type factor.
2.1. Governing Equations: The raw artificial-viscosity field is smoothed by a discrete averaging operator to produce a smooth continuous coefficient.
2.2. High-Order HDG Discretization: HDG couples element interiors through a single-valued mesh-skeleton trace, while curved high-order mappings represent physical boundaries and anisotropic boundary-layer meshes.
2.3. Implicit Temporal Discretization: DIRK advances the solution through diagonally implicit stages, and each stage retains the HDG algebraic structure with a conservative-state mass contribution.
2.5. Static Condensation: Static condensation reduces the coupled HDG system to a Schur-complement trace system, after which element corrections are recovered independently.

3. Massively Parallel HDG Solver

The solver exploits HDG hybridization to separate element-local dense algebra from globally coupled trace operations, combining additive-Schwarz-preconditioned GMRES with reduced-basis initialization for repeated condensed solves.

Solver architecture: Static condensation reduces the globally coupled HDG problem to a trace system, while conservative and gradient variables remain local and are recovered after trace correction.This decomposition underpins the solver’s parallel structure.
Solver architecture: HDG separates element-local dense linear algebra from global sparse communication, restricting interpartition communication to shared trace degrees of freedom.Local work includes dense factorizations, triangular solves, and matrix–matrix products associated with elemental Schur complements.
Linear solver: Restarted GMRES solves the condensed trace system, with a one-level additive Schwarz preconditioner exploiting local face-patch coupling and independent dense subdomain solves.The preconditioner includes no global coarse correction.
Linear solver: A reduced-basis approximation initializes GMRES for strongly correlated systems arising across Newton iterations, DIRK stages, and time steps.The basis is updated incrementally with recently accepted trace corrections and a fixed window of stored vectors.
Distributed GPU implementation: The implementation retains principal solution, residual, and operator data in device memory while managing intermediate data through a small scratchpad.The authors leave finer trade-offs among kernel fusion, data locality, and parallelism for future work.

4. Subsonic Transitional Flows over the Eppler 387 Airfoil

The Eppler 387 cases assess whether the HDG solver can resolve subsonic laminar–turbulent transition across four Reynolds numbers. Results show Reynolds-number-dependent upstream transition, shorter active regions and reattachment distances, and close agreement between predicted and measured pressure distributions.

4.1. Problem Definition: The validation cases examine subsonic transitional flow over Eppler 387 at four Reynolds numbers, focusing on separation, transition, reattachment, and aerodynamic response.The freestream Mach number is 0.09 and the angle of attack is 6°; diagnostics include pressure, skin friction, boundary-layer measures, fluctuations, and vortical structures.
4.2. Simulation Setup: The three-dimensional mesh contains 661,760 curved hexahedral elements from a 20,680-element C-grid extruded through 32 spanwise layers.Near-wall refinement targets the viscous sublayer, laminar boundary layer, separated shear layer, and near wake, with polynomial degree k = 2.
4.3. Transitional Flow Physics: As Reynolds number increases, three-dimensional vortices emerge farther upstream, develop over a shorter suction-side interval, and transition to turbulent reattachment more compactly.The highest-Reynolds-number case shows earlier three-dimensionality and a shorter streamwise distance to reattachment than the lowest-Reynolds-number case.
4.3. Transitional Flow Physics: Boundary-layer thickness measures and velocity fields show the largest and most extended separated-region deficit at the lowest Reynolds number.Displacement thickness grows through separation, while the low-speed layer contracts after reattachment; higher Reynolds numbers produce shorter recovery regions.
4.3. Transitional Flow Physics: Predicted pressure distributions closely reproduce McGhee et al.’s measurements across all four Reynolds numbers, including suction peak, pressure recovery, and suction-side plateau dependence.Skin friction independently indicates a broader near-zero or negative-Cf region at low Reynolds number and faster recovery to positive Cf at higher Reynolds numbers.
4.3. Transitional Flow Physics: The turbulent kinetic-energy and pressure-fluctuation regions shift upstream and contract with Reynolds number, remaining concentrated in the separated shear layer and near wake.The disturbance-amplitude footprint likewise becomes shorter and more closely attached to the suction-side separation region at higher Reynolds numbers.

5. Supersonic Taylor–Green Vortex

The supersonic Taylor–Green vortex benchmark tracks the development and decay of three-dimensional compressible turbulence, including transient shock-associated structures. Mesh refinement brings kinetic-energy and dissipation histories toward reference DNS behavior, with dilatational dissipation requiring especially fine resolution.

5.1. Problem Description: The benchmark uses a triply periodic cube with analytically specified initial conditions to isolate compressible turbulence development from wall-boundary effects.The flow begins as a smooth large-scale vortex field and develops finer vortices, compressive regions, shocks, and shock–turbulence interaction.
5.3. Compressible Turbulence: The solver examines kinetic energy, solenoidal dissipation, dilatational dissipation, and instantaneous flow structures as complementary diagnostics.Kinetic energy measures large-scale decay, solenoidal dissipation reflects vortical-cascade activity, and dilatational dissipation measures compressive shock-related activity.
5.3. Compressible Turbulence: The finest simulations evolve from organized large-scale vortices into a complex three-dimensional turbulent network with transient shocklet-like structures.At later times, pressure variations narrow and the strongest compressive events weaken as the turbulent field decays.
5.3. Compressible Turbulence: Mesh refinement produces monotonically convergent kinetic-energy histories, with the 963^3, 1123, and 1283 solutions closest to reference DNS during maximum dissipation.Coarse meshes lose energy too rapidly once the flow becomes strongly multiscale.
5.3. Compressible Turbulence: Solenoidal dissipation dominates overall, while dilatational dissipation peaks near 7 × 10^-4 and is roughly an order of magnitude smaller but more resolution-sensitive.The dilatational peak near t ≈ 6–7 coincides with the emergence of steep compressive structures, and fine meshes better reproduce reference behavior.
5.3. Compressible Turbulence: The benchmark shows faster convergence for kinetic energy than for dissipation measures as the flow transitions from organized vortices to compressible turbulence.The dilatational component remains particularly sensitive to the resolution of sharp compressive features.

6. Hypersonic Boundary-Layer Transition

The Mach-7 boundary-layer benchmark resolves curvature-driven instability, compressible disturbances, and nonlinear transition on a finely clustered high-order mesh. Increasing freestream turbulence advances breakdown, while weaker forcing preserves coherent streaks farther downstream and shifts transition-related wall responses downstream.

6.1. Problem Description: The hypersonic configuration combines thin boundary layers, strong compressibility, curved geometry, Görtler structures, acoustic disturbances, and nonlinear secondary breakdown.It therefore tests transition prediction in a strongly compressible wall-bounded flow across multiple instability mechanisms.
6.2. Simulation Setup: The production meshes contain more than 12 million quadratic elements, with near-wall clustering used to resolve hypersonic boundary-layer and thermal-layer dynamics.The simulations use k = 2 geometry and solution approximation, a three-stage third-order DIRK scheme, ∆t = 2 × 10^-4, and 9200 time steps.
6.2. Simulation Setup: The simulations start from a converged 2D laminar precursor mapped onto the 3D mesh, focusing subsequent evolution on imposed disturbances and centrifugal instability.The precursor removes the startup transient associated with forming the basic hypersonic boundary layer.
6.4. Boundary Layer Instability and Transition: Pressure and density fields progress downstream from smooth disturbances to increasingly complex fluctuations and strongly fluctuating transitional regions.The 0.5% turbulence-intensity case develops appreciable small-scale fluctuations farther upstream than the lower-intensity cases.
6.4. Boundary Layer Instability and Transition: Counterrotating streamwise vortices redistribute momentum and thermal energy, producing alternating enhanced and reduced wall heat-flux bands before nonlinear streak fragmentation.The 0.5% case loses streak coherence earliest, while the 0.2% and 0.1% cases retain organized structures farther downstream.
6.4. Boundary Layer Instability and Transition: Time-averaged heat-flux and skin-friction fields retain persistent streaks and shift the transitional region downstream as freestream turbulence intensity decreases.Lower-intensity forcing delays nonlinear breakdown while allowing curvature-driven streaks to develop coherently over a longer distance.
6.4. Boundary Layer Instability and Transition: Increasing freestream turbulence intensity accelerates transition without changing the qualitative curvature-driven transition pathway.Heat flux and skin friction provide complementary evidence for this shared behavior.
6.4. Boundary Layer Instability and Transition: The 0.5% case responds earlier upstream, whereas the 0.2% and 0.1% cases remain similar longer before their wall-coefficient increases and maxima occur farther downstream.The 0.1% heat-transfer response extends farther downstream than the 0.2% response.

7. Conclusions

The massively parallel GPU HDG solver combines high-order accuracy, shock robustness, implicit nonlinear solution, and scalable execution for compressible DNS. Across subsonic, supersonic, and hypersonic benchmarks, it agrees with experiments or published DNS while resolving transition, compressibility, shocks, and three-dimensional turbulence.

The solver combines high-order HDG, localized shock capturing, DIRK time integration, Newton–GMRES, additive Schwarz preconditioning, and reduced-basis acceleration.These components address accuracy, robustness, and scalability together for compressible DNS.
Eppler 387 simulations closely match experiments across four Reynolds numbers, capturing separation, instability, reattachment, and wake development.Agreement includes suction peaks, pressure recovery, and Reynolds-number-dependent pressure plateaus.
Supersonic Taylor–Green simulations converge toward reference DNS histories while resolving shocklet-like structures and three-dimensional turbulent dynamics.The mesh-sequence study examines kinetic-energy and dissipation histories.
Mach-7 simulations resolve Görtler-dominated transition through streak amplification, three-dimensional modulation, and nonlinear breakdown.Wall heat-flux and skin-friction distributions provide the clearest transition signatures.
Lowering freestream turbulence intensity from 0.5% to 0.2% and 0.1% shifts breakdown downstream while preserving the qualitative transition pathway.The dominant wall-response features shift downstream alongside the breakdown process.
The framework currently relies on carefully designed meshes and calorically perfect-gas equations, leaving adaptive resolution and real-gas thermochemistry as major extensions.Future work also targets adaptive time stepping, stronger preconditioners, reduced-order acceleration, and broader validation.