From h to p efficiently: Strategy selection for operator evaluation on hexahedral and tetrahedral elements

A spectral/hp element discretisation permits both geometric flexibility and beneficial convergence properties to be attained simultaneously. The choice of elemental polynomial order has a profound effect on the efficiency of different implementation strategies with their performance varying substantially for low and high order spectral/hp discretisations. We examine how careful selection of the strategy minimises computational cost across a range of polynomial orders in three dimensions and compare how different operators, and the choice of element shape, lead to different break-even points between the implementations. In three dimensions, higher expansion orders quickly lead to a large increase in the number of element-interior modes, particularly in hexahedral elements. For a typical boundary–interior modal decomposition, this can rapidly lead to a poor performance from a global approach, while a sum-factorisation technique, exploiting the tensor-product structure of elemental expansions, leads to better performance. Furthermore, increased memory requirements may cause an implementation to show poor runtime performance on a given system, even if the strict operation count is minimal, due to detrimental caching effects and other machine-dependent factors.