Source-linked AI summary
PyFR: An Open Source Framework for Solving Advection-Diffusion Type Problems on Streaming Architectures using the Flux Reconstruction Approach
Freddie D Witherden, Antony M Farrington, Peter E Vincent
TL;DR
High-fidelity CFD near complex geometries requires methods that combine high-order accuracy, unstructured-grid flexibility, and efficient execution on modern hardware. The paper presents PyFR, an open-source flux-reconstruction framework for mixed unstructured grids and streaming architectures, and demonstrates accuracy, GPU performance, and scalability. Its supported scope includes specific governing systems, element types, and interface assumptions.
Problem
Industrial CFD seeks high-fidelity scale-resolving simulations of transient compressible flows near complex geometries, while existing first- and second-order RANS software is not well suited to such simulations.
Method
PyFR is an open-source Python-based framework using flux reconstruction, a Mako-derived domain-specific language, and mixed unstructured grids targeted at CPUs and NVIDIA GPUs.
Results
PyFR shows spatial super accuracy for 2D Euler equations, expected order for Couette flow, and compute-bound kernels achieving 50%–90% of realisable peak FLOP/s.
Takeaways & Limitations
The demonstrated framework supports high-order simulations across mixed unstructured grids and streaming hardware, including qualitative unsteady 3D viscous-flow results on curved grids.
Takeaways & Limitations
The current release supports quadrilateral and triangular grids in two dimensions and hexahedral grids in three dimensions, while interface flux points must be co-located physically.
Abstract
from arXiv · showhide
High-order numerical methods for unstructured grids combine the superior accuracy of high-order spectral or finite difference methods with the geometric flexibility of low-order finite volume or finite element schemes. The Flux Reconstruction (FR) approach unifies various high-order schemes for unstructured grids within a single framework. Additionally, the FR approach exhibits a significant degree of element locality, and is thus able to run efficiently on modern streaming architectures, such as Graphical Processing Units (GPUs). The aforementioned properties of FR mean it offers a promising route to performing affordable, and hence industrially relevant, scale-resolving simulations of hitherto intractable unsteady flows within the vicinity of real-world engineering geometries. In this paper we present PyFR, an open-source Python based framework for solving advection-diffusion type problems on streaming architectures using the FR approach. The framework is designed to solve a range of governing systems on mixed unstructured grids containing various element types. It is also designed to target a range of hardware platforms via use of an in-built domain specific language based on the Mako templating engine. The current release of PyFR is able to solve the compressible Euler and Navier-Stokes equations on grids of quadrilateral and triangular elements in two dimensions, and hexahedral elements in three dimensions, targeting clusters of CPUs, and NVIDIA GPUs. Results are presented for various benchmark flow problems, single-node performance is discussed, and scalability of the code is demonstrated on up to 104 NVIDIA M2090 GPUs. The software is freely available under a 3-Clause New Style BSD license (see www.pyfr.org).
Program Description
PyFR is an open-source Python framework for solving advection-diffusion problems with high-order flux reconstruction on mixed unstructured grids and streaming hardware. Its release targets compressible flow equations across CPU and GPU systems.
- Framework: PyFR v0.1.0 is an open-source Python-based framework for solving compressible Euler and Navier-Stokes equations.The software is released under the New Style BSD license and uses Python, CUDA, and C.
- Hardware and parallelism: The framework targets clusters of CPUs and GPUs, with multi-GPU and multi-CPU awareness implemented through MPI and OpenMP.Its listed computing environment extends up to GPU clusters running recent Linux/UNIX systems.
- Numerical method: The solution method is a high-order flux reconstruction approach suitable for curved, mixed, unstructured grids.
- Implementation: PyFR uses symbolic manipulation and run-time code generation through a domain-specific language.The implementation relies on Python libraries including numpy, PyCUDA, mpi4py, SymPy, and Mako.
- Runtime: Many small problems can be solved on a recent workstation.
Nomenclature
The nomenclature defines indexing, element sets, polynomial expansions, coordinates, and operators used throughout the formulation. It establishes conventions for summation, dimensions, basis functions, and interface quantities.
- Index conventions: Repeated dummy indices on the right-hand side are summed, and all indices are zero-based.
- Functions: The notation defines common functions and matrix properties, including the Kronecker delta, determinant, and matrix dimensions.
- Elements: Element notation distinguishes element types, individual elements, the standard element, boundaries, and element counts.
- Expansions: Expansion notation specifies polynomial order, spatial dimensions, field variables, nodal basis polynomials, coordinates, and the transformed-to-physical mapping.
1 Introduction
The introduction motivates high-order CFD for transient compressible flows near complex geometries and explains why streaming architectures favor local, parallel computation. PyFR applies flux reconstruction and explicit time stepping to this setting.
- Motivation: Industrial CFD practitioners seek high-fidelity scale-resolving simulations of transient compressible flows near complex geometries.The cited examples include UAV spoilers, intake ducts, weapons bays, and complete configurations at Reynolds numbers 10^4–10^7 and Mach numbers M ∼0.1–1.0.
- Hardware context: Modern processors have increasingly outpaced memory bandwidth, raising the importance of memory locality for sustained computational throughput.The processor FLOPS-per-byte ratio increased from ∼0.2 to upwards of ∼4 between older and newer chips.
- Streaming computation: Streaming architectures divide computation into independent streams that can be distributed across cores and vector lanes.
- Numerical approach: Compute-intensive discretisations formulated for stream processing are well suited to current and likely future hardware platforms.The introduction identifies flux reconstruction combined with explicit time stepping as an archetypal example.
- PyFR objective: PyFR targets advection-diffusion problems on mixed unstructured grids across CPUs and NVIDIA GPUs using a Mako-derived domain-specific language.The current release supports quadrilateral and triangular elements in two dimensions and hexahedral elements in three dimensions.
2 Flux Reconstruction
Flux reconstruction maps fields onto standard elements, represents solutions and fluxes at nodal points, and couples neighboring elements through common interface values and fluxes. The resulting semi-discrete system is advanced in time after correction and divergence evaluation.
- Element representation: The formulation uses standard elements and mappings to transform physical coordinates and fields into a common computational space.Element types include quadrilaterals and triangles in two dimensions and several three-dimensional shapes.
- Point sets: Solution points define nodal basis functions, while flux points lie on element boundaries with associated outward normal vectors.The basis is constructed from a selected polynomial space using a generalized Vandermonde matrix.
- Interfaces: Flux points paired across an interface must share the same physical coordinates, ensuring geometric alignment between neighboring elements.Figure 2 illustrates this co-location for a triangle and quadrangle.
- Interface coupling: FR interpolates discontinuous solution values from solution points to flux points, then computes common interface solutions from neighboring values.On unstructured grids, an indirect lookup table identifies adjoining flux points, and common-solution functions may upwind or downwind values.
- Flux correction: Correction functions combine discontinuous polynomial derivatives with interface jumps to construct continuous-flux gradients and divergences.The common normal flux is defined to balance outward and inward fluxes across adjoining elements, preserving conservation.
- Time integration: The procedure produces a semi-discretised system of ordinary differential equations in time that can be advanced with schemes such as classical fourth-order Runge-Kutta.
3 Implementation
PyFR implements Flux Reconstruction through compact, platform-portable kernels and schedules computation to reuse memory and overlap communication with computation. Its Python-based domain-specific language generates CPU and NVIDIA GPU code while supporting distributed-memory execution.
- Overview: PyFR is a compact, efficient, platform-portable Python implementation of the Flux Reconstruction approach.
- Overview: The implementation separates constant operators, state matrices, matrix-multiply kernels, point-wise kernels, and distributed-memory scheduling.
- Operator and State Matrices: Symbolic manipulation with SymPy generates polynomial operators and derivatives at solution and flux points to keep the codebase compact.
- Operator and State Matrices: PyFR uses structure-of-arrays packing across platforms, with field-variable spacing equal to the number of solution points in an element.
- Matrix Multiplication: BLAS and vendor-optimized GEMM routines provide matrix multiplication, although PyFR's matrix sizes are not necessarily optimal for GEMM.
- Point-Wise Kernels: Mako-based runtime templates generate C/OpenMP and CUDA kernels, avoiding separate point-wise implementations for each target platform.
- Point-Wise Kernels: Point-wise kernels prescribe scalar operations while the templating engine handles vectorization and memory gathering; generated CUDA and OpenMP C variants target different hardware.
- Distributed Memory Parallelism and Scheduling: Persistent, non-blocking MPI communication and split interface kernels enable rank-local computation to overlap with ghost-state exchange.
4 Governing Systems
PyFR supports compressible Euler and Navier-Stokes governing systems within the Flux Reconstruction framework. The formulations define inviscid and viscous fluxes, interface treatments, and their two- or three-dimensional reductions.
- Overview: The current PyFR release solves compressible Euler equations for inviscid flow and compressible Navier-Stokes equations for viscous flow.
- Euler Equations: The three-dimensional Euler formulation uses conservative variables containing density, velocity, total energy, and pressure-related flux terms.
- Euler Equations: Euler interface fluxes are computed with approximate Riemann solvers, including the implementations described in PyFR's appendix.
- Compressible Navier-Stokes Equations: The Navier-Stokes flux is written as inviscid flux minus viscous flux, with viscous terms defined through stress-energy components and heat-transfer contributions.
- Compressible Navier-Stokes Equations: The thermal parameter Δ is defined as μC_p/Pr, where μ is dynamic viscosity and Pr is the Prandtl number.
- Compressible Navier-Stokes Equations: The LDG treatment computes common boundary solutions and augments the inviscid Riemann solver for viscous interface fluxes, with β controlling upwinding or downwinding.
- Dimensional Reduction: The three-dimensional equations reduce to two dimensions by deleting the fourth variable rows and third flux columns, yielding two-dimensional velocity vectors.
5 Validation
PyFR was validated on Euler vortex and compressible Navier–Stokes benchmarks using structured, mixed unstructured, extruded hexahedral, and fully unstructured hexahedral meshes. The tests recovered expected accuracy, including super accuracy for the Euler vortex, while a turbulent cylinder wake was simulated at scale.
- Euler Equations: Euler Vortex Super Accuracy: The Euler vortex test used four structured quadrilateral meshes with spacings h = 1/3, 2/7, 1/4, and 2/9.The vortex had strength S = 13.5, free-stream Mach number M = 0.4, and radius R = 1.5.
- Euler Equations: Euler Vortex Super Accuracy: FR schemes DG, SD, and HU exhibited super accuracy beyond the expected ℘+1 order for the Euler vortex.The error initially reflects projection error before the time-dependent spatial-operator error dominates.
- Compressible Navier-Stokes Equations: Couette Flow: The Couette-flow solver obtained the expected order of accuracy on four two-dimensional mixed meshes containing triangles and quadrilaterals.Accuracy was assessed using experimental L2 energy errors.
- Compressible Navier-Stokes Equations: Couette Flow: The expected order of accuracy was obtained on extruded hexahedral meshes and fully unstructured hexahedral meshes.The fully unstructured test used meshes with 96, 536, and 1004 elements, although those results had higher standard errors.
6 Single Node Performance
Single-node evaluation on an NVIDIA M2090 characterized PyFR’s matrix and point-wise kernels across polynomial orders. Matrix-multiplication performance remained high, while weak scalability was evaluated through simulations reaching 10^4 GPUs.
- Single GPU Performance: The evaluation used a fully periodic mesh containing 50,176 hexahedral elements and Navier–Stokes polynomial orders ℘ = 2, 3, and 4.Kernel categories were DGEMM, direct-memory point-wise kernels, and indirectly accessed point-wise kernels.
- Single GPU Performance: 305.4 GFLOP/s was achieved by the ℘ = 4 matrix-multiplication kernels, exceeding 75% of the realisable peak.The absolute performance regressed slightly at ℘ = 4, increasing the fraction of wall-clock time spent in these kernels.
- Single GPU Performance: Matrix-multiplication time increased with polynomial order, whereas point-wise-kernel arithmetic performance remained invariant with order.Changing order increases the number of processed points without changing the point-wise operation itself.
7 Scalability
PyFR scalability is evaluated for three-dimensional Navier–Stokes simulations on NVIDIA M2090 GPUs, covering both strong and weak scaling. Strong scaling is near perfect through eight GPUs, while weak scaling reaches 104 GPUs for a 3.8 × 10^9-degree-of-freedom case.
- Hardware and setup: The Emerald cluster combines NVIDIA M2090 GPUs across 84 nodes connected by QDR InfiniBand.The cluster includes nodes with either three or eight M2090 GPUs.
- Hardware and setup: A ℘= 3 Navier–Stokes case uses 114 688 structured hexahedral elements and approximately 4720 MiB of memory.This working set loads an M2090 to approximately 90% of its user-available memory when ECC is enabled.
- Weak scalability: 3.8 × 10^9 degrees of freedom are simulated at N = 104 in the weak-scaling study, with a working set of approximately 485 GiB.The domain is enlarged along the y-axis to preserve the decomposition strategy as GPU count increases.
- Strong scalability: Strong-scaling speedup is near perfect through eight GPUs, then degrades; at N = 32, the speedup is 26.At 32 GPUs, each GPU is loaded below 3%, explaining the reduced scaling efficiency.
8 Conclusions
The conclusions present PyFR as an accurate and performant framework for advection-diffusion problems on streaming architectures. The reported demonstrations span accuracy, GPU kernel performance, and scalability for Navier–Stokes simulations.
- Framework: PyFR is an open-source Python framework for solving advection-diffusion type problems on streaming architectures.Its targeting methodology supports multiple hardware platforms.
- Accuracy and flow demonstrations: PyFR shows spatial super accuracy for 2D Euler equations and expected accuracy for Couette flow across 2D and 3D grids.The conclusions also report qualitative results for unsteady 3D viscous flows on curved grids.
- Performance: Compute-bound kernels reach 50%–90% of realisable peak FLOP/s, while bandwidth-bound point-wise kernels exceed 89% of realisable peak bandwidth.These performance figures are reported for an NVIDIA M2090 GPU in three dimensions.
- Scalability: Strong scalability is demonstrated to 32 NVIDIA M2090s and weak scalability to 104 NVIDIA M2090s for the 3D Navier–Stokes equations.Both scalability results concern the same GPU architecture.
A Matrix Representation
The matrix representation recasts most Flux Reconstruction operations as matrix–matrix multiplications. Constant operator matrices and packed state matrices provide a compact algebraic form for the governing equations and their implementation.
- Matrix formulation: Most operations in an FR step can be expressed as matrix–matrix multiplications.This form exposes the element-local algebra used by the implementation.
- Matrix formulation: The reformulation introduces a constant operator matrix A and state matrices B and C, with scalar coefficients c1,2.The operator matrix is fixed, whereas B and C contain solution-dependent state data.
- Data packing: Field variables for each element may be packed along matrix rows in multiple ways, creating an implementation degree of freedom.The selected packing determines how subsequent matrix operations traverse element data.
- Operator construction: The same matrix procedure is applied to the formulations corresponding to Equations 6, 9, and 11.The construction includes context-dependent strides and a block-diagonal structure for M5.
B Approximate Riemann Solvers
The framework represents interfaces through common solution and flux prescriptions, including boundary-specific ghost states and approximate Riemann solvers. Boundary types are parameterized according to the physical condition being imposed.
- Interface fluxes: At an interface, uL and uR are the discontinuous solution states, with the normal vector associated with the first state defining the flux orientation.The right inviscid flux is evaluated from uR, and interface inviscid fluxes are prescribed separately.
- Approximate Riemann solvers: The local Lax–Friedrichs, or Rusanov, solver prescribes inviscid numerical interface fluxes using an estimate s of the maximum wave speed.The wave-speed estimate enters the approximate Riemann flux construction.
- Boundary conditions: Boundary interfaces use one flux point and boundary-type functions to prescribe common solutions and normal fluxes.The functions receive the solution, gradient, and unit normal at the boundary flux point.
- Boundary conditions: Boundary conditions may instead provide ghost solution states and gradients through a boundary transformation B(b).Different ghost states can be supplied for the common solution and normal-flux functions, and the ghost gradient may depend on uL.
- Boundary-condition parameters: Supersonic inflow uses free-stream density, velocity, and pressure, whereas subsonic outflow uses free-stream pressure.The no-slip isothermal wall depends on wall temperature and wall velocity, usually with vw = 0.