# Convergence problem solving steady Navier-Stokes with multiphenicsx

Hello everyone.
I am trying to find a solver for the steady Navier-Stokes equations that I need to solve in a complex 2D domain with high velocity values at the inlet.
Therefore, I found this useful tutorial: tutorial_navier_stokes to which I am referring, and to do a first tentative I simply modified the constant inlet velocity increasing its value from 1 to 100, but by doing so the non-linear solver does not converge, stopping after two iterations and not calculating correctly the velocity field.
The part of the tutorial I modified is simply:

nu = 0.01

def u_in_eval(x: np.typing.NDArray[np.float64]) -> np.typing.NDArray[  # type: ignore[no-any-unimported]
petsc4py.PETSc.ScalarType]:
"""Return the flat velocity profile at the inlet."""
values = np.zeros((2, x.shape[1]))
values[0, :] = 100.0
return values

def u_wall_eval(x: np.typing.NDArray[np.float64]) -> np.typing.NDArray[  # type: ignore[no-any-unimported]
petsc4py.PETSc.ScalarType]:
"""Return the zero velocity at the wall."""
return np.zeros((2, x.shape[1]))

And when I run the solver as :

up_m = run_monolithic()
(u_m, p_m) = (up_m.sub(0).collapse(), up_m.sub(1).collapse())

I get:

0 542.3530856245577
1 538.968125845104

Then the solver stops, and when plotting the velocity field I can notice the values are wrong and the solver did not converge.
Why this solver does not work for high values of the inlet? Maybe the Reynolds becomes too high and this solver does not work for high Reynolds number? How can I solve this issue?
For completeness, I copy and paste a MWE code to get this result:

import typing
import basix.ufl
import dolfinx.fem
import dolfinx.fem.petsc
import dolfinx.io
import gmsh
import mpi4py.MPI
import numpy as np
import numpy.typing
import petsc4py.PETSc
import ufl
import viskex
import multiphenicsx.fem
import multiphenicsx.fem.petsc
#constitutive parameters
nu = 0.01
def u_in_eval(x: np.typing.NDArray[np.float64]) -> np.typing.NDArray[  # type: ignore[no-any-unimported]
petsc4py.PETSc.ScalarType]:
"""Return the flat velocity profile at the inlet."""
values = np.zeros((2, x.shape[1]))
values[0, :] = 100.0
return values

def u_wall_eval(x: np.typing.NDArray[np.float64]) -> np.typing.NDArray[  # type: ignore[no-any-unimported]
petsc4py.PETSc.ScalarType]:
"""Return the zero velocity at the wall."""
return np.zeros((2, x.shape[1]))

# geometrical parameters
pre_step_length = 4.
after_step_length = 14.
pre_step_height = 3.
after_step_height = 5.
mesh_size = 1. / 5.

#Mesh

gmsh.initialize()
p0 = gmsh.model.geo.addPoint(0.0, after_step_height - pre_step_height, 0.0, mesh_size)
p1 = gmsh.model.geo.addPoint(pre_step_length, after_step_height - pre_step_height, 0.0, mesh_size)
p2 = gmsh.model.geo.addPoint(pre_step_length, 0.0, 0.0, mesh_size)
p3 = gmsh.model.geo.addPoint(pre_step_length + after_step_length, 0.0, 0.0, mesh_size)
p4 = gmsh.model.geo.addPoint(pre_step_length + after_step_length, after_step_height, 0.0, mesh_size)
p5 = gmsh.model.geo.addPoint(0.0, after_step_height, 0.0, mesh_size)
l0 = gmsh.model.geo.addLine(p0, p1)
l1 = gmsh.model.geo.addLine(p1, p2)
l2 = gmsh.model.geo.addLine(p2, p3)
l3 = gmsh.model.geo.addLine(p3, p4)
l4 = gmsh.model.geo.addLine(p4, p5)
l5 = gmsh.model.geo.addLine(p5, p0)
line_loop = gmsh.model.geo.addCurveLoop([l0, l1, l2, l3, l4, l5])
gmsh.model.geo.synchronize()
gmsh.model.addPhysicalGroup(1, [l0, l1, l2, l4], 2)
gmsh.model.mesh.generate(2)
mesh, subdomains, boundaries = dolfinx.io.gmshio.model_to_mesh(
gmsh.model, comm=mpi4py.MPI.COMM_WORLD, rank=0, gdim=2)
gmsh.finalize()
# Create connectivities required by the rest of the code
mesh.topology.create_connectivity(mesh.topology.dim - 1, mesh.topology.dim)
boundaries_1 = boundaries.indices[boundaries.values == 1]
boundaries_2 = boundaries.indices[boundaries.values == 2]
#Function Spaces
V_element = basix.ufl.element("Lagrange", mesh.basix_cell(), 2, shape=(mesh.geometry.dim, ))
Q_element = basix.ufl.element("Lagrange", mesh.basix_cell(), 1)
#Standard formulation using a mixed function space
def run_monolithic() -> dolfinx.fem.Function:
"""Run standard formulation using a mixed function space."""
# Function spaces
W_element = basix.ufl.mixed_element([V_element, Q_element])
W = dolfinx.fem.functionspace(mesh, W_element)

# Test and trial functions: monolithic
vq = ufl.TestFunction(W)
(v, q) = ufl.split(vq)
dup = ufl.TrialFunction(W)
up = dolfinx.fem.Function(W)
(u, p) = ufl.split(up)

# Variational forms
F = (nu * ufl.inner(ufl.grad(u), ufl.grad(v)) * ufl.dx
+ ufl.inner(ufl.grad(u) * u, v) * ufl.dx
- ufl.inner(p, ufl.div(v)) * ufl.dx
+ ufl.inner(ufl.div(u), q) * ufl.dx)
J = ufl.derivative(F, up, dup)

# Boundary conditions
W0 = W.sub(0)
V, _ = W0.collapse()
u_in = dolfinx.fem.Function(V)
u_in.interpolate(u_in_eval)
u_wall = dolfinx.fem.Function(V)
u_wall.interpolate(u_wall_eval)
bdofs_V_1 = dolfinx.fem.locate_dofs_topological((W0, V), mesh.topology.dim - 1, boundaries_1)
bdofs_V_2 = dolfinx.fem.locate_dofs_topological((W0, V), mesh.topology.dim - 1, boundaries_2)
inlet_bc = dolfinx.fem.dirichletbc(u_in, bdofs_V_1, W0)
wall_bc = dolfinx.fem.dirichletbc(u_wall, bdofs_V_2, W0)
bc = [inlet_bc, wall_bc]

# Class for interfacing with SNES
class NavierStokesProblem:
"""Define a nonlinear problem, interfacing with SNES."""

def __init__(  # type: ignore[no-any-unimported]
self, F: ufl.Form, J: ufl.Form, solution: dolfinx.fem.Function,
bcs: list[dolfinx.fem.DirichletBC], P: typing.Optional[ufl.Form] = None
) -> None:
self._F = dolfinx.fem.form(F)
self._J = dolfinx.fem.form(J)
self._obj_vec = dolfinx.fem.petsc.create_vector(self._F)
self._solution = solution
self._bcs = bcs
self._P = P

def create_snes_solution(self) -> petsc4py.PETSc.Vec:  # type: ignore[no-any-unimported]
"""
Create a petsc4py.PETSc.Vec to be passed to petsc4py.PETSc.SNES.solve.

The returned vector will be initialized with the initial guess provided in `self._solution`.
"""
x = self._solution.vector.copy()
with x.localForm() as _x, self._solution.vector.localForm() as _solution:
_x[:] = _solution
return x

def update_solution(self, x: petsc4py.PETSc.Vec) -> None:  # type: ignore[no-any-unimported]
"""Update `self._solution` with data in `x`."""
with x.localForm() as _x, self._solution.vector.localForm() as _solution:
_solution[:] = _x

def obj(  # type: ignore[no-any-unimported]
self, snes: petsc4py.PETSc.SNES, x: petsc4py.PETSc.Vec
) -> np.float64:
"""Compute the norm of the residual."""
self.F(snes, x, self._obj_vec)
return self._obj_vec.norm()  # type: ignore[no-any-return]

def F(  # type: ignore[no-any-unimported]
self, snes: petsc4py.PETSc.SNES, x: petsc4py.PETSc.Vec, F_vec: petsc4py.PETSc.Vec
) -> None:
"""Assemble the residual."""
self.update_solution(x)
with F_vec.localForm() as F_vec_local:
F_vec_local.set(0.0)
dolfinx.fem.petsc.assemble_vector(F_vec, self._F)
dolfinx.fem.apply_lifting(F_vec, [self._J], [self._bcs], x0=[x], scale=-1.0)
dolfinx.fem.set_bc(F_vec, self._bcs, x, -1.0)

def J(  # type: ignore[no-any-unimported]
self, snes: petsc4py.PETSc.SNES, x: petsc4py.PETSc.Vec, J_mat: petsc4py.PETSc.Mat,
P_mat: petsc4py.PETSc.Mat
) -> None:
"""Assemble the jacobian."""
J_mat.zeroEntries()
dolfinx.fem.petsc.assemble_matrix(  # type: ignore[misc]
J_mat, self._J, self._bcs, diagonal=1.0)  # type: ignore[arg-type]
J_mat.assemble()
if self._P is not None:
P_mat.zeroEntries()
dolfinx.fem.petsc.assemble_matrix(  # type: ignore[misc]
P_mat, self._P, self._bcs, diagonal=1.0)  # type: ignore[arg-type]
P_mat.assemble()

# Create problem
problem = NavierStokesProblem(F, J, up, bc)
F_vec = dolfinx.fem.petsc.create_vector(problem._F)
J_mat = dolfinx.fem.petsc.create_matrix(problem._J)

# Solve
snes = petsc4py.PETSc.SNES().create(mesh.comm)
snes.setTolerances(max_it=20)
snes.getKSP().setType("preonly")
snes.getKSP().getPC().setType("lu")
snes.getKSP().getPC().setFactorSolverType("mumps")
snes.setObjective(problem.obj)
snes.setFunction(problem.F, F_vec)
snes.setJacobian(problem.J, J=J_mat, P=None)
snes.setMonitor(lambda _, it, residual: print(it, residual))
up_copy = problem.create_snes_solution()
snes.solve(None, up_copy)
problem.update_solution(up_copy)  # TODO can this be safely removed?
up_copy.destroy()
snes.destroy()
return up
up_m = run_monolithic()
(u_m, p_m) = (up_m.sub(0).collapse(), up_m.sub(1).collapse())
viskex.dolfinx.plot_vector_field(u_m, "u")

Two suggestions:

1. Please compute the Reynolds number with the values of viscosity and inlet velocity as in the original tutorial. I can’t remember it off the top of my head (and don’t want to recompute it myself right now), but I would bet it was something between 50 and 100.
2. Then, compute the Reynolds number in your modified case. If it is 100 times larger than 50~100, i.e. 5e3~1e4, than you’ll likely need some stabilization and/or a better initial guess.

I cannot find the original tutorial where the Reynolds number it’s computed, so I do not know which parameter was used for the geometrical reference, but yes it should be 100 times larger than 50-100.
However, even by setting the inlet velocity equal to 50 ( so Reynolds maximum of the order 5e3), I still get these residuals:

0 271.17654281227885
1 269.92777904571636

And also by setting the inlet velocity equal to 20 ( Reynolds max around 2e3) I cannot reach convergence even with 100 iterations.
Also, it is not very clear to me how to set a different initial guess in this code, if you could maybe help me understanding it.
However, I would need a solver for the Navier Stokes that can handle also Reynolds number of the order of 1e3, do you maybe have any suggestions about it?
Thank you very much.

I cannot find the original tutorial where the Reynolds number it’s computed

You literally linked it in your first post. I can see inlet velocity magnitude equal to 1, viscosity equal to 0.01 and inlet height (say we use this as reference length) equal to 3, which makes Re = 300.

And also by setting the inlet velocity equal to 20

Which would have Re = 6e3: that is a “large” value of Re, and in general you cannot expect a plain FE method to handle it. Possible suggestions:
0) you need some physical understanding of the fluid dynamics problem. Is the flow going to be steady or unsteady? Is the flow going to be laminar or “turbulent” (in quotes, because to properly speak about turbulence you would need a 3D flow)

1. you need some stabilization, for instance SUPG. Based on 0), in the case of the backward facing step, you may also need some backflow stabilization if the recirculation region is wider than the horizontal length behind the step
2. based on 0), you need to understand if it makes sense to keep using a stationary simulation, or if you need to go for an unsteady one. Honestly, I doubt that the underlying physical phenomenon is actually steady with Re of the order of O(1e3).