This work discusses the development of a three-dimensional, high-order, compressible viscous flow solver for mixed unstructured grids that can run on multiple GPUs. The solver utilizes a range of so-called Vincent-Castonguay-Jameson-Huynh (VCJH) flux reconstruction schemes in both tensor-product and simplex elements. Such schemes are linearly stable for all orders of accuracy and encompass several well known high-order methods as special cases. Because of the high arithmetic intensity associated with VCJH schemes and their element-local nature, they are well suited for GPUs. The single-GPU solver developed in this work achieves speed-ups of up to 45x relative to a serial computation on a current generation CPU. Additionally, the multi-GPU solver scales well, and when running on 32 GPUs achieves a sustained performance of 2.8 Teraops (double precision) for 6th-order ac- curate simulations with tetrahedral elements. In this paper, the techniques used to achieve this level of performance are discussed and a performance analysis is presented. To the authors' knowledge, the aforementioned flow solver is the first high-order, three-dimensional, compressible Navier-Stokes solver for mixed unstructured grids that can run on multiple GPUs.