Indeed with higher velocity the solver diverges so that SUPG/PSPG is needed.
SUPG/PSPG won’t necessarily help with nonlinear convergence, but it may help with the accuracy of the resulting converged solution. For difficult nonlinear problems, you might look into using the PETSc SNES solver with a backtracking line search. This can be accessed through FEniCS’s Python API as shown in this discussion.
Then the solver becomes extremely slow, diverges and I get this message below.
By default, FEniCS will try to estimate an appropriate polynomial degree to integrate exactly with Gauss quadrature. However, for complicated variational forms (such as those obtained from stabilized methods), this estimate is usually overkill. For complicated problems, it’s usually a good idea to manually set the degree of numerical quadrature, e.g.,
dx = dx(metadata={"quadrature_degree":2*k})
(and similar for ds and dS, if applicable) which is what the warning message is suggesting to do.
Also: Another type of stabilization that can add robustness for incompressible Navier–Stokes is “least squares on the incompressibility constraint” (LSIC; called “grad-div” in some references), which is also defined in my notes. It essentially penalizes the deviations from \nabla\cdot\mathbf{u} = 0 that can occur at individual points when incompressibility is only enforced weakly with respect to some finite-dimensional test space.