Cache Blocking for Flux Reconstruction: Extension to Navier-Stokes Equations and Anti-aliasing

TL;DR

In this article, cache blocking is implemented for the Navier Stokes equations with anti-aliasing support on mixed grids in PyFR for CPUs to eliminate unnecessary data movements between kernels at the main memory level.

Abstract

In this article, cache blocking is implemented for the Navier Stokes equations with anti-aliasing support on mixed grids in PyFR for CPUs. In particular, cache blocking is used as an alternative to kernel fusion to eliminate unnecessary data movements between kernels at the main memory level. Specifically, kernels that exchange data are grouped together, and these groups are then executed on small sub-regions of the domain that fit in per-core private data cache. Additionally, cache blocking is also used to efficiently implement a tensor product factorisation of the interpolation operators associated with anti-aliasing. By using cache blocking, the intermediate results between application of the sparse factors are stored in per-core private data cache, and a significant amount of data movement from main memory is avoided. In order to assess the performance gains a theoretical model is developed, and the implementation is benchmarked using a compressible 3D Taylor-Green vortex test case on both hexahedral and prismatic grids, with third- and forth-order solution polynomials. The expected performance gains based on the theoretical model range from 1.99 to 2.62, and the speedups obtained in practice range from 1.67 to 3.67 compared to PyFR v1.11.0.

Figures (10)

• Figure 1: Mapping function from physical space to reference space for a quadrilateral element.
• Figure 2: Solution points and flux points for a quadrilateral and a triangular elements.
• Figure 3: Kernel prerequisites map for Navier-Stokes Solver including optional Flux AA kernels.
• Figure 4: Current kernel execution order. Arrows going over green line indicate main memory access.
• Figure 5: Solution points for prism elements with polynomial orders $p=2$ and $p=3$. Point distribution is based on William-Shun at the triangular cross sections and based on Gauss-Legendre along the z-direction. A direct interpolation between $p=2$ and $p=3$ prisms solution points results in a dense operator matrix as any point in $p=3$ depend on all the points in $p=2$ as indicated with orange points.