DSLHyPE: Bilingual MLIR Compiler for PDE Kernels

• DSLHyPE is a bilingual Python-embedded DSL and compiler that separates numerical schemes from native C/C++ PDE physics and generates CPU or GPU kernels through MLIR.
• It supports patch-based explicit finite-volume and finite-difference methods, including stencils, halo regions, Riemann solvers, source terms, adaptive meshes, and nonlinear hyperbolic PDE systems.
• Experiments with CCZ4 and Euler solvers show faster generated code than basic C++ baselines, reduced numerical-loop code size, and major memory-locality benefits from flattening nested data layouts.

DSLHyPE is a bilingual domain-specific language and compiler for specifying numerical compute kernels within ExaHyPE, the Exascale Hyperbolic PDE Engine. Its defining separation is between the numerical scheme, expressed in a Python-embedded DSL, and PDE physics, retained in existing native C or C++ implementations or generated from SymPy. DSLHyPE lowers the numerical description and native physics into MLIR, applies compiler transformations, and generates CPU or GPU implementations. The system is designed for hyperbolic PDE applications involving finite volumes, finite differences, discontinuous Galerkin-related machinery, explicit time stepping, adaptive meshes, nonlinear systems, custom Riemann solvers, admissibility checks, and computationally intensive source or flux terms (Stokes et al., 18 Aug 2026).

1. Design objectives and computational scope

DSLHyPE addresses the entanglement of numerical algorithms and physical models in high-performance PDE software. Numerical researchers often need to change stencils, reconstruction procedures, flux differencing, source-term treatment, Runge–Kutta stages, or dissipation mechanisms, while existing physical implementations may contain substantial tested C or C++ code. Rewriting both the numerical scheme and the physics for every target architecture creates duplicated CPU and GPU implementations and complicates maintenance.

The system separates two layers:

1. Numerical scheme and update structure: expressed in DSLHyPE and comprising stencils, data ranges, halo accesses, temporary arrays, reconstruction, flux differencing, source-term application, and time-integration stages.
2. PDE physics: implemented in native C or C++, including conservative fluxes, source terms, non-conservative products, eigenvalue calculations, Riemann solvers, admissibility checks, nonlinear constitutive relations, and boundary-condition calculations.

The underlying PDE class is represented by

$\partial_t Q+\nabla F(Q)+B(Q)\cdot\nabla Q=S(Q),$

where $Q\in\mathbb{R}^{N+M}$ contains $N$ evolved quantities and $M$ fixed material or model parameters, $F(Q)$ is the conservative flux, $B(Q)\cdot\nabla Q$ denotes non-conservative products, and $S(Q)$ is a source term.

DSLHyPE is not a replacement for ExaHyPE’s generic infrastructure or for existing PDE implementations. It provides a compiler-visible numerical representation around them. ExaHyPE supplies finite-volume, finite-difference, discontinuous Galerkin-related, explicit time-stepping, adaptive-mesh, and CPU/GPU capabilities, while DSLHyPE specifies and compiles the cell- or patch-level numerical kernel.

2. Kernel abstraction and bilingual programming model

For a single cell, ExaHyPE represents a compute kernel as a mapping

$$\mathcal{K}_{\Delta t}: \mathbb{R}^{\hat p^d(N+M)} \longrightarrow \mathbb{R}^{p^dN},$$

where $d$ is the spatial dimension, $p$ is the number of interior points per dimension, $Q\in\mathbb{R}^{N+M}$0 is the halo width, and $Q\in\mathbb{R}^{N+M}$1 is the extent including halo points. The input contains evolved and material quantities over the halo-extended patch, whereas the output contains the evolved variables over the interior region.

For a two-dimensional patch, the logical input extent is approximately

$Q\in\mathbb{R}^{N+M}$2

while the output covers the $Q\in\mathbb{R}^{N+M}$3 evolved variables over a $Q\in\mathbb{R}^{N+M}$4 interior. ExaHyPE generally batches multiple cells:

$Q\in\mathbb{R}^{N+M}$5

where $Q\in\mathbb{R}^{N+M}$6 is the number of cells processed together. In the reported experiments, Peano’s adaptive mesh traversal can assemble batches of up to $Q\in\mathbb{R}^{N+M}$7 cells.

A kernel is viewed as a composition of numerical stages,

$Q\in\mathbb{R}^{N+M}$8

A finite-volume implementation may evaluate source terms, add scaled source contributions, reconstruct states at faces, evaluate fluxes or Riemann solvers, apply flux differences, add dissipation, and produce the updated cell. DSLHyPE represents this staged structure through assignments to DataBlock objects.

The language is bilingual:

• DSLHyPE/Python specifies arrays, ranges, halo regions, shifted accesses, temporary data, and the sequence of numerical operations.
• C/C++ supplies the physical functions used by the numerical scheme.
• MLIR provides the common intermediate representation into which both sides are lowered.

This arrangement permits existing physics implementations to remain in their native language while allowing the numerical scheme to be optimized independently.

3. LLM and data abstractions

DSLHyPE is embedded in Python. Kernel functions remain syntactically valid Python and are parsed using Python’s ast machinery, but their bodies use DSLHyPE abstractions rather than ordinary Python execution.

A representative declaration is:

$M$1

Quoted type names identify compiled native structures. A runtime value can be bound into the DSL using:

$M$2

Symbols such as N, M, p, h, dt, and dx may be global or statically configured parameters.

DataBlock

The primary abstraction is a DataBlock, representing an array with explicit logical bounds and a source or destination expression. An input patch can be declared as:

$M$3

Its dimensions represent the component index, two spatial indices, and the batch-cell index. Negative lower bounds express halo coordinates directly, allowing stencil operations to use logical rather than manually offset indices.

An output region can be declared as:

$M$4

The paper identifies three forms of data block:

• Existing-array aliases, which refer to arrays owned by the surrounding application.
• Function-evaluation blocks, whose values are produced by applying a native function over a declared range.
• Computed temporary blocks, formed through stencil operations, pointwise expressions, or global matrix applications.

A function-evaluation block can represent a flux calculation:

$M$5

The interface identifies the native routine, its arguments, and the expected output shape. Native functions may compute fluxes, source terms, maximum wave speeds, Riemann fluxes, or other physics quantities.

Slicing and shifted indexing

Python-like slicing denotes logical regions:

$M$6

Shifted slices express neighboring stencil accesses. A representative update is:

$M$7

The shifted regions represent adjacent face or stencil values. DSLHyPE lowers these operations into indexed loads and stores in generated loops.

Assignments to data blocks imply iteration over the relevant ranges. The user does not explicitly write CPU threads, GPU blocks, or SIMD lanes. Intermediate data blocks expose stage dependencies, giving the compiler visibility into producers, consumers, read regions, and write regions.

4. Compilation and MLIR integration

DSLHyPE’s compilation pipeline combines Python AST processing, MLIR generation, memory-layout transformation, native-code translation, and target-specific lowering:

$Q\in\mathbb{R}^{N+M}$9

The principal stages are the following.

Python frontend

Python’s ast module is used to identify DataBlock declarations, ranges, slices, assignments, arithmetic expressions, native-function calls, and runtime bindings. The system has both a pretty printer that produces baseline C++ and an MLIR-generating compiler.

DSL MLIR generation

The DSL compiler lowers ranged array operations into standard MLIR dialects, including:

• arith for arithmetic;
• scf for structured control flow and loops;
• memref for allocation and memory references;
• GPU and LLVM-related dialects during later lowering.

At this stage, numerical stages appear as loop-oriented operations over declared data ranges.

Native physics translation

C and C++ physics code is translated into MLIR with Polygeist. This avoids implementing an independent C/C++ parser and makes the native routines available to the same optimization pipeline as the DSL-generated code.

The two representations must be integrated before final lowering. The custom pass

$M$8

combines DSL-generated and Polygeist-generated IR. Early integration permits cross-function optimization, inlining, whole-kernel vectorization, and GPU device lowering. If the physics remains a separately compiled host function, these transformations cannot inspect the complete numerical kernel.

Nested-memory linearization

ExaHyPE commonly stores cell data using an array-of-structs organization, with batches represented effectively as arrays of arrays of structs. A direct MLIR representation can resemble nested memory references such as

$M$9

or conceptually a double**. Such indirection can impair vectorization, cache behavior, GPU address calculation, and dependence analysis.

The custom pass

$F(Q)$0

linearizes this representation by gathering input data into a contiguous region, executing the flattened kernel, and scattering results back. The pass transforms nested references into a flat memref<?x?xf64> representation and allocates temporary storage as required.

GPU memory-space transformation

The pass

$F(Q)$1

rewrites memory operations for device execution. It extends the gather and scatter stages to allocate device memory, copy data to the GPU, execute device kernels, copy results back, and release or manage temporary storage.

Target lowering

CPU lowering proceeds through LLVM-oriented representations to x86 code. GPU lowering uses MLIR’s GPU dialect and produces host launch code, device kernels, and embedded PTX/device binaries. Runtime callbacks such as

$F(Q)$2

provide the GPU runtime interface; on the Grace-Hopper system, these callbacks are implemented using CUDA.

The generated DSLHyPE IR is intended to remain independent of a specific GPU architecture. A different target, such as an AMD GPU, would require different device-lowering passes and runtime callbacks rather than a different numerical DSL representation.

5. Numerical schemes, memory layout, and optimization

DSLHyPE is demonstrated with explicit finite-volume and finite-difference schemes. The numerical computation uses Cartesian patch data, halo regions, stencils, pointwise algebra, and time-level updates. The compiler exposes the regular loop structure that would otherwise remain hidden behind runtime abstractions.

The system supports numerical operations involving conservative fluxes, source terms, non-conservative products, Riemann solvers, reconstructed face states, dissipation, directional derivatives, and matrix or tensor-like operations. Physics functions may contain nonlinear expressions, branches, local variables, and mathematical-library calls.

Specialization

Static or configured quantities can become compile-time values, including:

• patch dimensions;
• halo widths;
• field and component dimensions;
• stencil ranges;
• loop bounds;
• spatial offsets;
• time-step and grid-spacing parameters.

This specialization allows constant propagation, loop simplification, inlining, and target-specific lowering.

Memory organization

ExaHyPE’s native representation is commonly array-of-structs, with each cell’s data contiguous. A batch of patches is effectively an array of arrays of structs. DSLHyPE’s linearization pass converts this representation into flat contiguous storage around the core kernel.

The flattened representation reduces nested pointer indirection. In the reported CCZ4 configuration with patch size $N$0 and $N$1 patches:

• unflattened MLIR transferred approximately $N$2 of the main-memory data volume of native C++;
• flattening reduced this to approximately $N$3;
• unflattened MLIR increased cycles spent in D-TLB load misses by approximately $N$4;
• its L2-cache memory-bound metric increased by almost $N$5;
• with flattening, the L2 difference fell to approximately $N$6.

These measurements indicate that memory representation, rather than arithmetic generation alone, is a major determinant of generated-kernel performance.

Optimization opportunities

The current implementation preserves the DSL’s sequence of stages relatively conservatively. The paper identifies several optimization opportunities:

• Kernel fusion: combining stages to reduce temporary-memory traffic and launch overhead.
• AoS-to-SoA conversion: transforming data layout for SIMD or GPU execution while keeping the physics code unchanged.
• Temporary reuse: allocating shared scratch storage based on data lifetimes rather than creating a separate temporary for every stage.
• Loop fusion and permutation: improving cache and TLB locality.
• Tensor and linear-algebra lowering: recognizing linear stencil sections and mapping them to linear-algebra dialects or BLAS-like routines.
• Automatic cost models: selecting fusion, batching, layout, and target strategies rather than relying on a manually selected pass sequence.

The current integration uses pure kernel offloading in ExaHyPE. Cell data are not retained persistently on the GPU, so batches may incur gathering, allocation, host-to-device transfer, kernel launch, device-to-host transfer, and scattering. These costs are especially important for small kernels and schemes with large halo regions.

6. Evaluation, applications, and limitations

DSLHyPE is evaluated with two PDE systems:

• CCZ4 gravitational-wave solver: 59 nonlinear PDE variables, a non-conservative formulation, source terms, and finite-volume and fourth-order finite-difference variants.
• Euler matter-evolution solver: five equations, a conservative flux, and finite-volume and fourth-order finite-difference variants.

The experiments use three-dimensional patches of sizes

$N$7

and compare four configurations: Euler finite volume, Euler fourth-order finite difference, CCZ4 finite volume, and CCZ4 fourth-order finite difference.

The hardware and software environment comprises an Intel Xeon Platinum 8480+ Sapphire Rapids processor, a Grace-Hopper H200 GPU with 141 GB of available memory, LLVM/MLIR 21.1.8, and LIKWID 5.5.1. Measurements use a standalone kernel driver and average ten runs. The comparisons include C++ generated by a simple pretty printer, C++ with OpenMP, DSLHyPE-generated MLIR with or without memory flattening, and GPU MLIR.

CPU performance

Across the reported CPU configurations, MLIR-generated implementations are faster than the vanilla C++ implementations. Larger patches generally provide higher throughput. Increasing the batch size $N$8 initially improves performance for small patches, but excessive batching eventually increases runtime. For large patches, performance is less sensitive to $N$9. CCZ4 is less sensitive to surrounding loop structure because its nonlinear physics functions dominate execution time.

The paper does not provide absolute seconds-per-cell tables in the supplied text; the relevant values are presented graphically. It does report that generated MLIR code improves over the basic generated C++ baseline and that memory flattening substantially reduces the overhead associated with nested memory representations.

GPU performance

The GPU experiments compare DSLHyPE GPU MLIR with C++ generated using OpenMP target offload. DSLHyPE performs uniformly better in the tested configurations. Compute time grows approximately linearly with workload after an initial small-$M$0 burn-in phase. Higher-order schemes can achieve better pure compute utilization after GPU saturation because their larger arithmetic workload provides more work per launch.

Total runtime is substantially larger than pure device compute time because it includes launch, gathering, allocation, and host-device transfers. For fourth-order finite differences, larger halo layers increase data movement. The paper reports that kernel launches, data reorganization, and data migration can approximately double the higher-order scheme’s runtime, particularly for the MLIR variant.

Code size

The reported core numerical-loop line counts are:

Scheme	DSLHyPE	Handwritten CPU	Handwritten GPU
FV	85	660	1124
FD4	73	752	Not implemented

These counts exclude user-defined Euler and CCZ4 physics code. They show that one numerical specification can replace separately maintained CPU and GPU numerical-loop implementations.

Principal limitations

DSLHyPE’s current limitations include:

• incomplete GPU mathematical-library support, with functions such as fmax and fmin potentially preventing successful lowering;
• difficulty with extensive dynamic stack allocation during MLIR inlining;
• manually selected optimization passes rather than automatic fusion, layout, batching, and target decisions;
• GPU offloading overhead from nonpersistent device data;
• temporary-memory allocation and reuse costs;
• limited support for explicitly specified reductions, synchronization, user-defined parallelism, and general control flow;
• no complete standalone grammar in the supplied description;
• no formal type system reported for the DSL;
• no detailed syntax for boundary-condition declarations or automatic differentiation;
• no quantitative compilation-time breakdown;
• no comprehensive convergence, error-norm, or independent numerical-validation study in the supplied evaluation.

The demonstrated numerical scope is primarily explicit finite-volume and finite-difference computation on patch-based data. Broader support for adaptive numerical strategies, more sophisticated tensor transformations, improved memory reuse, and automatic optimization remains future work.

DSLHyPE’s significance lies in making the numerical algorithm compiler-visible without requiring researchers to rewrite established PDE physics. Its bilingual model preserves native C or C++ implementations at the source level, combines them with Python-expressed numerical schemes during compilation, and lowers the unified representation through MLIR to CPU or GPU code. The reported results demonstrate reduced numerical-loop code volume, x86 and H200 generation, improved performance over basic generated C++, and successful integration with gravitational-wave and matter-evolution solvers. The principal remaining challenge is to replace manually engineered lowering choices with compiler-driven decisions for fusion, data layout, temporary reuse, batching, scheduling, and persistent accelerator execution (Stokes et al., 18 Aug 2026).