A matrix-free stabilized solver for the incompressible Navier-Stokes equations

We present an efficient solver for the incompressible Navier-Stokes equations implemented in a matrix-free fashion. It uses a higher-order continuous Galerkin finite element method for the space discretization and leverages a stabilized formulation that includes both the SUPG and PSPG terms. We solve the non-linear problem in a fully coupled way, using a Newton-Krylov method, which is preconditioned by a monolithic geometric multigrid solver. To evaluate its efficiency in terms of time to solution and scalability on modern high-performance computers, we use a manufactured solution, a steady flow around a sphere with Reynolds number Re = 150 and the Taylor–Green vortex benchmark at Re = 1 600. The results indicate that the solver is robust and scales for both steady-state and transient problems. We compare the matrix-free solver to a matrix-based version and show it exhibits lower memory requirements, better scalability, and significant speedups (10–100× for higher-order elements). Moreover, we demonstrate that a matrix-free implementation is highly efficient when using higher-order elements, which provide higher accuracy at a lower number of degrees of freedom for complex steady problems. To the best of our knowledge, this work is the first that uses a matrix-free monolithic geometric multigrid preconditioner to solve the stabilized Navier-Stokes equations. All implementations are available via the open-source software Lethe.