Architecture-aware h-to-p optimisation: spectral/hp element operators for mixed-element meshes
Comput. Phys. Commun., vol. 330, pp. 110467 (2027)
@article{xing-2027,
title = {Architecture-aware $h$-to-$p$ optimisation: spectral/$hp$ element operators for mixed-element meshes},
author = {Xing, J. Y. and Xia, B. and Renner, D. and Cantwell, C. D. and Moxey, D. and Kirby, R. M. and Sherwin, S. J.},
journal = cpc,
volume = {330},
pages = {110467},
year = {2027},
doi = {10.1016/j.cpc.2026.110467},
url = {https://www.sciencedirect.com/science/article/pii/S0010465526004480},
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.}
}
Spectral element codes are usually tuned for hexahedral elements, but a mesh of a real geometry also contains prisms, pyramids and tetrahedra. This paper works out how the standard finite element operators should be implemented for each of those shapes, finding that the best strategy depends on the element shape, the polynomial order and the machine. It also sets out a way of evaluating operators such as the stiffness matrix that makes more use of the collocation properties of the expansion. On a GPU the Helmholtz operator on tetrahedra runs at most 2.5 times slower than on hexahedra, despite needing six times as many floating-point operations.
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.