Mixed-Precision in High-Order Methods: the Impact of Floating-Point Precision on the ADER-DG Algorithm

TL;DR

The paper presents a mixed-precision ADER-DG implementation in ExaHyPE, enabling per-kernel precision control to study accuracy across five hyperbolic PDE scenarios. It demonstrates that full high-order convergence generally requires fp64 precision for all kernels, while fp32 supports moderate orders and fp16/bf16 frequently fail or produce unstable results. Mixed-precision configurations can mitigate some failures, but gains are highly problem- and kernel-dependent, with predictor often driving accuracy outcomes. The work lays out a template-based code-generation approach to implement and test mixed precision in a production-ready PDE engine, providing guidance on kernel selection and precision trade-offs for accuracy-critical simulations.

Abstract

We present a mixed-precision implementation of the high-order discontinuous Galerkin method with ADER time stepping (ADER-DG) for solving hyperbolic systems of partial differential equations (PDEs) in the hyperbolic PDE engine ExaHyPE. The implementation provides a simple API extension for specifying the numerical precision for individual kernels, and thus allows for testing the effect of low and mixed precision on the accuracy of the solution. To showcase this, we study the impact of precision on the overall convergence order and actual accuracy of the method as achieved for four common hyperbolic PDE systems and five relevant scenarios that feature an analytic solution. For all scenarios, we also assess how sensitive each kernel of the ADER-DG algorithm is to using double, single or even half precision. This addresses the question where thoughtful adoption of mixed precision can mitigate hurtful effects of low precision on the overall simulation.

Figures (12)

• Figure 1: Example how to generate C++ template specialisations in all required precisions using the Jinja2 template engine.
• Figure 2: Different declaration of the user function flux(), depending on whether the user has specified multiple precisions (right) or to perform all computations in uniform double precision (left). Note that flux() is called by the predictor and the corrector.
• Figure 3: Initial density $\rho(x,y) = 0.02 \left(1+e^{-50(x^2+y^2)}\right)$ of the Gaussian bell scenario for the Euler equations.
• Figure 4: Initial density $\rho(x,y) = \rho_{\infty} + \delta \rho$ of the static isentropic vortex scenario for the Euler equations.
• Figure 5: L2 error for the initial solution of the elastic planar-wave scenario (Sec. \ref{['subsubseq:planar_waves']}) with different polynomial orders, cell sizes and precisions.