High-order numerical methods for unstructured grids combine the superior accuracy of high-order spectral or finite difference methods with the geometrical 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 many-core hardware platforms, such as graphics 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 talk we will present PyFR an open-source, Python-based framework for solving advection-diffusion type problems using the FR approach. The framework is designed to target a range of hardware platforms via use of a domain specific language. With this PyFR is able to solve the compressible Euler and Navier-Stokes equations on grids of quadrilateral and triangular elements in two dimensions, and hexahedral, tetrahedral, prismatic, and pyramidal elements in three dimensions, targeting clusters of multi-core CPUs, NVIDIA GPUs, AMD GPUs, Intel Xeon Phis, and heterogeneous mixtures thereof. Results will be presented for various benchmark and "real-world" flow problems, and scalability/performance of PyFR will be demonstrated on clusters with thousands of NVIDIA GPUs. Throughout the talk the importance of algorithm-software-hardware co-design, in the context of next-generation computational fluid dynamics, will be highlighted.