There are several points to consider:
-
First, the Bubnov–Galerkin formulation is inherently unstable at a high element Reynolds number (i.e., the Reynolds number taking the element diameter as a length scale), so it will become impractical to use at a high global Reynolds number (which would force a very small element size for stability). This affects the accuracy of the solution even if you can coax the nonlinear solver into converging, but it can also make the equation system harder to solve. The remedy for this is some form of stabilization, e.g., SUPG/PSPG, also considered a type of variational multiscale (VMS) method. (See Chapter 6 of my course notes.)
-
Second, the algebraic residual is difficult to converge, even for stabilized formulations! The usual checklist of things to try is:
- Use a line search scheme, which can be accessed by using a PETSc SNES through a
NonlinearVariationalSolver(as demonstrated in FEniCS here). - Use a continuation method, e.g., solve once with a low Reynolds number, use that solution as the initial guess for a slightly higher Reynolds number, etc.
- If nothing else works, solve an unsteady problem and try to drive it to a steady state. (However, this may not always work well for Navier–Stokes, because if turbulence develops it may never settle into a steady solution.)
- Use a line search scheme, which can be accessed by using a PETSc SNES through a