Architecture-aware $h$-to-$p$ optimisation: spectral/$hp$ element operators for mixed-element meshes

Abstract: We extend earlier international efforts to optimise hexahedral-based spectral element methods on GPUs and vectorised CPUs to mixed element meshes additionally involving prismatic, pyramidic, and tetrahedral shapes using tensorial expansions. We demonstrate that common finite element operators (such as the mass and Helmholtz matrices) benefit from alternative implementation strategies depending on the element shape, choice of polynomial order, and system architecture in order to achieve optimal performance. In addition, we introduce a new approach/interpretation to efficiently evaluate more complex operations involving inner products with the derivative of the expansions as part of the integrand such as the stiffness matrix. This approach seeks to maximise operations using the collocation properties of the nodal tensorial expansion associated with classical quadrature rules. Our GPU performance tests demonstrate that the throughput of the Helmholtz operator on tetrahedral elements is at most 2.5 times slower than on hexahedral elements, despite tetrahedra having a factor of six greater floating-point operations.

• The paper presents an architecture-aware framework that extends h-to-p optimisation for spectral/hp element operators to unstructured, mixed-element meshes.
• It employs tensorial expansions and tailored implementation strategies (StdMat, SumFac, SumFacTOP) to maximize performance on both CPUs and GPUs.
• Empirical benchmarks demonstrate that operator reformulation, including collocation-based methods, significantly enhances throughput across varied element types and polynomial degrees.

Architecture-aware $h$-to-$p$ Optimisation for Spectral/$hp$ Element Operators on Mixed-Element Meshes

Introduction

The development of high-order finite element methods (FEM), specifically spectral/$hp$ element techniques, for PDE solvers on unstructured, mixed-element meshes is driven by the need for higher accuracy and geometric flexibility in large-scale simulations across disciplines such as computational fluid dynamics. The increasing arithmetic intensity of $hp$-FEM aligns well with modern SIMD architectures—namely GPUs and vectorised CPUs—by improving FLOP-per-byte ratios and thus exploiting hardware design trends in high-performance computing. Efficient operator applications on such architectures demand not only mathematical but also architecture-aware algorithmic adaptation to element type, polynomial degree ($P$), and memory access. This paper extends established optimisations for hexahedral spectral element operators to mixed-element (including prismatic, pyramidic, and tetrahedral) settings via tensorial expansions, with a focus on optimally mapping finite element operators to heterogeneous computational platforms.

Mathematical and Algorithmic Framework

High-order mixed-element spectral/$hp$ FEM discretise PDEs by expanding the solution in a tailored set of basis functions (modal or nodal, polynomial up to order $P$) over elements of various shapes. The canonical example is the Helmholtz equation, discretised into an elemental matrix problem comprising both mass and stiffness (Helmholtz) operators. Hexahedral and quadrilateral elements allow true tensor-product expansions, enabling classic sum-factorisation techniques. Simplex-type elements (triangles, tetrahedra, pyramids, prisms) require collapsed tensor-product expansions, typically implemented through Duffy transformations to the reference element, with resulting integration and basis evaluation complexities.

Operator application boils down to a sequence of transformations between coefficient and physical (quadrature) space (codified by basis evaluation matrices), differentiations, and inner products. For non-affine transformation elements, increased quadrature order is required for spectral accuracy. Critically, the alternative modal-to-nodal transformation viewpoint for operator evaluation, wherein the action in quadrature space is sandwiched between basis transforms (pre- and post-multiplied by the modal-to-nodal and its transpose), allows maximal use of collocation properties and potentially reduced floating-point operation counts on SIMD hardware.

For hybrid elements on unstructured meshes, tensorisation on the collapsed coordinate system is carefully addressed, especially for triangles and their 3D analogs, allowing the exploitation of sum-factorisation even for non-tensorial element types.

Architecture-aware Implementation Strategies

The paper systematises implementation strategies in three principal dimensions: memory space models, execution space models, and per-operator algorithmic choices, all realised in Nektar++. The memory model distinguishes HostSpace (CPU/vector) and DeviceSpace (GPU) allocations. Execution spaces, following a design inspired by Kokkos, are mapped hierarchically (serial, SIMD vectorised, device accelerators), with appropriate unification of backend APIs.

Block-based data structures coalesce same-type elements (shape, $P$, quadrature, mapping) for batch execution and memory efficiency. Operator calls are dispatched per block, allowing per-block algorithmic specialisation (e.g., matrix-based, sum-factorisation, threaded, or collocation-based versions).

Three principal operator implementations are considered:

• StdMat (Standard Matrix): Matrix-based, small-block GEMM, optimal for low $P$ or memory-bandwidth-bound use, especially in the mass operator and for elements with minimal tensor structure.
• SumFac (Sum-Factorisation): Tensor contraction using 1D basis operations, highly advantageous for high $p$0 elements with strong tensor-product structure.
• SumFacTOP (Sum-Factorisation Threaded On output Points): Exploits device-level thread groups to parallelise over output points, providing scalable fine-grained parallelism and memory reuse for shared-memory accelerators (notably GPUs).

Figures 3 and 4 visualise the differing parallelism/memory access in SumFac and SumFacTOP, highlighting device-specific memory layout and thread assignment.

Figure 2: Diagram of SumFac GPU implementation, with contiguous memory access by threads within work groups, each operating on an element group.

Figure 4: Diagram of SumFacTOP GPU implementation, where thread groups distribute work across elements with parallel output point computation.

Operator formulations for the Helmholtz matrix are presented in both non-collocated and collocated forms; the latter exploits basis collocation to improve computational efficiency on SIMD hardware.

Performance Analysis on CPU and GPU Architectures

Systematic benchmarking is performed on contemporary heterogeneous platforms: NVIDIA GH200 Grace Hopper Superchip (GPU) and Intel Xeon 6526Y (AVX-512 CPU). Throughput is reported as degrees of freedom processed per second for mass and Helmholtz operator applications, with $p$1 and element shape varied across tests.

Key observations include:

• GPU: SumFacTOP delivers superior performance on the Helmholtz operator for all shapes and $p$2, exploiting device-level shared memory and parallelism; for the mass operator, StdMat and SumFacTOP are comparable at higher $p$3, but SumFac excels for simplex-type elements at low $p$4 due to efficient register usage. Notably, the throughput reduction from hexahedral to tetrahedral elements for Helmholtz is at most $p$5 (with a $p$6 increase in underlying FLOPs), underscoring the efficacy of the architecture-aware strategies.

Figure 1: Throughput for the mass operator on GH200 (GPU) and Xeon (CPU) platforms across element types and $p$7, showing best configurations for each regime.

Figure 3: Throughput for the Helmholtz operator on GH200 (GPU) and Xeon (CPU) across element types and $p$8.

• CPU: StdMat dominates at low $p$9; SumFac overtakes as $hp$0 increases, especially for non-hexahedral elements where memory and cache behaviour become limiting. The characteristic performance drop-off at large degrees of freedom evidences cache saturation effects.
• Affine vs. Deformed Elements: Affine (regular) elements benefit from reduced geometric data movement, leading to higher throughput especially for the Helmholtz operator, an effect pronounced on CPUs with limited cache. This is quantified in Figures 7 and 8.

Figure 5: Xeon throughput for the mass operator (SumFac) with regular vs. deformed geometry.

Figure 6: Xeon throughput for the Helmholtz operator (SumFac) with regular vs. deformed geometry.

• Collocation-based Optimisation: Rewriting operations to fully exploit collocation in the basis functions (collocated approach) yields $hp$1–$hp$2 performance improvement for hexahedra and $hp$350% for tetrahedra compared to the classical approach, as shown in Figure 9.

Figure 7: Xeon throughput for the Helmholtz operator with collocated vs. non-collocated SumFac implementations, illustrating significant gains of the collocated scheme.

Data Structures and Software Design

The field/block abstraction aligns memory layout and operator batching for SIMD/execution model requirements. This is illustrated in Figure 10, encapsulating Field and Block association for contiguous memory allocation and efficient parallel access.

Figure 8: Schematic of Field and Block data structure model for efficient contiguous device access in mixed-element computations.

Discussion and Implications

This work provides robust empirical evidence that implementation strategies for high-order FEM operators must be adapted to both the polynomial degree and element type, as well as the execution hardware. The relative performance advantages of StdMat, SumFac, and SumFacTOP are highly context-dependent—a single implementation cannot serve all regimes efficiently. Architecture-aware operator dispatch and block-level data organisation are shown to be essential for performance portability, especially on complex mixed-element, industrially relevant meshes.

A notable highlight is the contradictory throughput behaviour: despite increased computational complexity, the performance gap between hexahedral and tetrahedral elements under optimal strategies is surprisingly narrow.

Future implications include:

• The demonstrated performance portability mechanisms are applicable to evolving hardware architectures, including many-core CPUs and next-generation accelerators.
• The explicit use of operator reformulation (especially collocation-based approaches) has substantial impact and could guide the design of future matrix-free and high-order PDE solver frameworks.
• Block-wise, field-structured data models will underpin scalable implementations for exascale platforms, with direct consequences for load balancing and adaptive mesh refinement in industrial contexts.

Conclusion

The study establishes that architecture-aware algorithmic specialisation is mandatory for high performance with spectral/$hp$4 element operators on unstructured, mixed-element meshes. Efficiently leveraging tensor-product basis functions, block execution, and collocation properties across heterogeneous architectures delivers substantial throughput gains, ensuring that increased algorithmic complexity in hybrid-element discretisations does not preclude advantageous hardware utilisation. The presented strategies set a reference point for future spectral/$hp$5 FEM code design, highlighting the inseparability of numerical algorithm choices and system architecture for scalable, portable simulation software.