When dealing with the incompressible Navier-Stokes equations, the solution for the pressure is usually the most time-consuming kernel. In this work, we present recent developments to enhance the pressure convergence and scalability of our new in-house, open-source, high-order CFD code, SOD2D, designed from scratch to run fully on GPUs. Its main objective is to perform scale-resolving large eddy simulation on complex geometries that require correct energy conservation, which is ensured by a low-dissipation spectral element method (SEM) discretisation. An implicit-explicit BDF-EXT3 third-order operator splitting approach is used to solve the velocity-pressure coupling. The fractional splitting allows us to avoid pressure stability issues when equal-order interpolation of pressure and velocity is used. It requires solving a Laplace equation for the pressure at each time step, which is typically the most time-consuming kernel. Moreover, un- less multigrid is used, the number of iterations grows as the mesh is refined, inhibiting correct weak scalability. In our previous code, Alya, we have obtained excellent algorithmic scalability up to 3000 million DoF using the open-source Algebraic Multigrid library PSCToolkit. When moving to SOD2D, the high-order discretisation requires a matrix-free implementation for computational efficiency, which precludes the direct use of an AMG library. We are therefore following the approach suggested in [1], where a Laplacian matrix is assembled for a linear mesh with exactly the same degrees of freedom as the high-order mesh. Algebraic multigrid applied to the linear Laplacian matrix serves as an excellent preconditioner for the high-order Laplacian which is solved with a conjugate gradient solver using a matrix-free implementation. For the moment we have tested the implementation on meshes up to 150 million DoF using up to 32 Nvidia H100 GPUs on the MareNostrum V supercomputer. Convergence has improved substantially compared to our previous implementation that relied on a diagonally preconditioned conjugate gradient solver for the pressure. For the presentation we will include results on the weak scalability to further highlight the advantages of using multigrid.