In this talk I will discuss the technical and design challenges of developing high-order finite element flow solvers which can efficiently use graphics processing units. I will first describe the Galerkin-Boltzmann model for nearly incompressible flows and detail how to increase overall efficiency by using multirate/semi-analytic time integration and absorbing boundary layers. I will then outline the components of the libParanumal incompressible flow solver that included GPU accelerated preconditioners and asses the computationally intensive parts. I will detail how we used empirical roofline models to optimize the performance of the most time-consuming kernels by tuning the fine-grain parallelism, memory utilization and maximizing achieved device memory bandwidth. I will present some computational results obtained in two-and-three dimensional flows.