Can Explicit Subgrid Models Enhance Implicit LES Simulations? A GPU-Oriented High-Order-Solver Perspective

TL;DR

This work addresses how explicit subgrid-scale (SGS) models interact with high-order discontinuous Galerkin methods in turbulence simulations on GPUs. By systematically varying split forms, Riemann solvers, and the Vreman SGS model for the Taylor–Green vortex at Re=1600 and in the inviscid limit, it shows that implicit DG dissipation suffices in well-resolved regimes, while explicit SGS is beneficial in under-resolved cases but must be carefully tuned to avoid degrading intermediate scales. A key finding is the absence of a single configuration that delivers optimal accuracy across laminar, transitional, and turbulent phases, underscoring the need for adaptive, scale-aware stabilization (e.g., SVV or data-driven SGS tuning). The results provide practical guidance for selecting numerical strategies in high-order turbulence simulations and point toward adaptive approaches for next-generation GPU-based solvers. This has broad implications for efficiently leveraging very high polynomial orders in DG methods while maintaining spectral fidelity and stability.

Abstract

High-order Discontinuous Galerkin (DG) methods offer excellent accuracy for turbulent flow simulations, especially when implemented on GPU-oriented architectures that favor very high polynomial orders. On modern GPUs, high-order polynomial evaluations cost roughly the same as low-order ones, provided the DG degrees of freedom fit within device memory. However, even with high-order discretizations, simulations at high Reynolds numbers still require some level of under-resolution, leaving them sensitive to numerical dissipation and aliasing effects. Here, we investigate the interplay between intrinsic DG dissipation mechanisms (implicit dissipation) -- in particular split forms and Riemann solvers -- and explicit subgrid-scale models in Large Eddy Simulations (LES). Using the three-dimensional Taylor--Green vortex at $Re = 1600$ and an inviscid case ($Re \to \infty$), we evaluate kinetic energy dissipation, spectral accuracy, and numerical stability. Our results show that when stability for under-resolved turbulence is ensured through split-forms (energy- or entropy-stable) schemes, subgrid-scale (SGS) LES models are not strictly necessary. At moderate Reynolds numbers, when the spatial resolution is sufficient to capture the relevant turbulence scales (i.e., in well-resolved LES), adding SGS models does not improve accuracy because the wavenumber range where they act overlaps with the inherent numerical dissipation of the DG scheme. In contrast, when the resolution is insufficient, as is typical at high Reynolds numbers, explicit subgrid-scale models complement the numerical dissipation and enhance accuracy by removing the excess energy that numerical fluxes alone cannot dissipate. These findings provide practical guidance for choosing numerical strategies in high-order turbulence simulations.

Figures (10)

• Figure 1: GPU efficiency of HORSES3D on a single GPU. Efficiency increases with $P$, roughly doubling from $P=3$ to $P=7$, due to higher arithmetic intensity and reduced memory traffic.
• Figure 2: Kinetic energy dissipation rate for standard (non-split) formulations with Gauss nodes. Both central and Roe-based fluxes eventually diverge due to aliasing-driven instabilities, while the split-form scheme introduces stabilizing dissipation.
• Figure 3: Effect of central split forms on dissipation and spectral behavior. The standard and Morinishi schemes are shown only in Fig. \ref{['fig:part1_kinenrate']}, as they became unstable before $t/t_c=9$. Results for all schemes are shown in both figures, and the curves overlap, indicating only minimal differences between the methods. Overall, central split forms exhibit excess dissipation during transition and insufficient dissipation once turbulence is fully developed.
• Figure 4: Effect of the Vreman SGS model on the Chandrasekhar split form. The model does not improve the accuracy of the transitional dynamics, while in the fully turbulent regime larger constants recover the expected spectral decay at high wavenumbers. Increasing the constant shifts model-induced dissipation toward lower wavenumbers, overdamping intermediate scales, while smaller constants under-dissipate the highest wavenumbers.
• Figure 5: Comparison between the baseline central scheme and the same scheme augmented with a Vreman SGS model of different strengths for the TGV at Reynolds 1600. Each subfigure shows the energy spectra (top) and the spectral difference $\Delta E(k)$ (bottom). The figures highlight how the SGS model constant controls both the magnitude of the added dissipation and the wavenumber at which it becomes active.