Stokes problem with only Dirichlet boundary condition

Good morning,
I am trying to solve the Stokes problem with full Dirichlet boundary condition. In this case the problem is ill posed, because the pressure is defined up to a constant. To resolve the issue I am trying to add the null space. In the code below, what I am trying to do is to attach to the matrix a vector composed of zeros for the dofs relative to the velocity and ones for the dofs relative to the pressure. The problem is that somehow the compiler says that this is not the right basis.
Am I missing something?

import numpy as np
import dolfinx
import ufl
from petsc4py import PETSc
from mpi4py import MPI


N = 2
mesh = dolfinx.mesh.create_unit_cube(MPI.COMM_WORLD, N, N, N,
                                     cell_type=dolfinx.mesh.CellType.hexahedron)
P2 = ufl.VectorElement("Lagrange", mesh.ufl_cell(), 2)
P1 = ufl.FiniteElement("Lagrange", mesh.ufl_cell(), 1)
TH = ufl.MixedElement([P2, P1])
W = dolfinx.fem.FunctionSpace(mesh, TH)

def u_ex(x): 
    return np.vstack((4*(x[1]-1)*x[1]*(x[0]-1)*(x[0]-1)*x[0]*x[0]*(x[1]-0.5), -4*(x[0]-0.5)*x[1]*x[1]*x[0]*(x[1]-1)*(x[1]-1)*(x[0]-1), np.zeros_like(x[0])))
def p_ex(x):
    return x[0]+x[1]+x[2]-3./2

V, V_to_W = W.sub(0).collapse()
Q, Q_to_W = W.sub(1).collapse()

### Define the boundary condition -> which in this case are only of Dirichlet type 
uD = dolfinx.fem.Function(V)
uD.interpolate(u_ex)
# Create facet to cell connectivity required to determine boundary facets
tdim = mesh.topology.dim
fdim = tdim - 1
mesh.topology.create_connectivity(fdim, tdim)
boundary_facets = dolfinx.mesh.exterior_facet_indices(mesh.topology)
boundary_dofs = dolfinx.fem.locate_dofs_topological(V, fdim, boundary_facets)
bcs = [dolfinx.fem.dirichletbc(uD, boundary_dofs)]

(u, p) = ufl.TrialFunctions(W)
(v, q) = ufl.TestFunctions(W)

x = ufl.SpatialCoordinate(mesh)
f = ufl.as_vector([1-24*(-1 + x[0])*(x[0]-1)*x[0]*x[0]*(-0.5 + x[1])- 48*(0.166667 - x[0] + x[0]*x[0])*(-1 + x[1])*(-0.5 + x[1])*x[1],
                1+24*(-0.5 + x[0])*(-1 + x[1])*(x[1]-1)*x[1]*x[1] + 8*(-1 + x[0])*(-0.5 + x[0])*x[0]*(1 - 6*x[1] + 6*x[1]*x[1]),
                1.0+0*x[0]])

dx = ufl.Measure("dx", domain=mesh)
a = ufl.inner(ufl.grad(u), ufl.grad(v)) * dx + \
    ufl.div(v)*p*dx + q*ufl.div(u)*dx
L = ufl.inner(f, v) * dx

a_f = dolfinx.fem.form(a)
L_f = dolfinx.fem.form(L)
A = dolfinx.fem.petsc.assemble_matrix(a_f, bcs=bcs)
A.assemble()

# Tp print the matrix
from scipy.sparse import csr_matrix
assert isinstance(A, PETSc.Mat)
ai, aj, av = A.getValuesCSR()
Asp = csr_matrix((av, aj, ai))
print(Asp)

b = dolfinx.fem.petsc.assemble_vector(L_f)
dolfinx.fem.petsc.apply_lifting(b, [a_f], [bcs])
b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
dolfinx.fem.petsc.set_bc(b, bcs)
b.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)

# Create vector that spans the null space
# To see dof: https://fenicsproject.discourse.group/t/how-to-obtain-dofs-of-functionspace-in-dolfinx/7090
num_dofs_global_u = V.dofmap.index_map.size_local * V.dofmap.index_map_bs
num_dofs_global = W.dofmap.index_map.size_local * W.dofmap.index_map_bs
null_vec = np.zeros(num_dofs_global)
print(num_dofs_global, ", ", null_vec.size)
null_vec[num_dofs_global_u:] = 1.0
print(null_vec)

def array2petsc4py(g):
    Xpt = PETSc.Vec().create()
    print("*")
    Xpt.setSizes(g.shape[0])
    print("*")
    # Allocate memory for the PETSc vector
    Xpt.setUp()
    print("*")
    # Copy the values from the NumPy array to the PETSc vector
    Xpt.array = g
    Xpt.assemble()
    return Xpt

null_vec_petsc = array2petsc4py(null_vec)
null_vec_petsc.normalize()
nullspace = PETSc.NullSpace().create(vectors=[null_vec_petsc])
assert nullspace.test(A)

A.setNearNullSpace(nullspace)

Here is a minimal example highlighting the two issues of your code

  1. The boundary condition is wrong
  2. The nullspace creation is overly complicated
import dolfinx
import numpy as np
import ufl
from mpi4py import MPI
from petsc4py import PETSc

N = 10
mesh = dolfinx.mesh.create_unit_cube(MPI.COMM_WORLD, N, N, N,
                                     cell_type=dolfinx.mesh.CellType.hexahedron)

P2 = ufl.VectorElement("Lagrange", mesh.ufl_cell(), 2)
P1 = ufl.FiniteElement("Lagrange", mesh.ufl_cell(), 1)
TH = ufl.MixedElement([P2, P1])
W = dolfinx.fem.FunctionSpace(mesh, TH)

def u_ex(x): 
    return np.vstack((4*(x[1]-1)*x[1]*(x[0]-1)*(x[0]-1)*x[0]*x[0]*(x[1]-0.5), -4*(x[0]-0.5)*x[1]*x[1]*x[0]*(x[1]-1)*(x[1]-1)*(x[0]-1), np.zeros_like(x[0])))
def p_ex(x):
    return x[0]+x[1]+x[2]-3./2

V, V_to_W = W.sub(0).collapse()
Q, Q_to_W = W.sub(1).collapse()

### Define the boundary condition -> which in this case are only of Dirichlet type 
uD = dolfinx.fem.Function(V)
uD.interpolate(u_ex)
# Create facet to cell connectivity required to determine boundary facets
tdim = mesh.topology.dim
fdim = tdim - 1
mesh.topology.create_connectivity(fdim, tdim)
boundary_facets = dolfinx.mesh.exterior_facet_indices(mesh.topology)
boundary_dofs = dolfinx.fem.locate_dofs_topological((W.sub(0),V), fdim, boundary_facets)
bcs = [dolfinx.fem.dirichletbc(uD, boundary_dofs, W.sub(0))]

(u, p) = ufl.TrialFunctions(W)
(v, q) = ufl.TestFunctions(W)

dx = ufl.Measure("dx", domain=mesh)
a = ufl.inner(ufl.grad(u), ufl.grad(v)) * dx + \
    ufl.div(v)*p*dx + q*ufl.div(u)*dx

a_f = dolfinx.fem.form(a)
A = dolfinx.fem.petsc.assemble_matrix(a_f, bcs=bcs)
A.assemble()
ns_vec = dolfinx.fem.Function(W)
ns_vec.x.array[Q_to_W] = 1
dolfinx.la.orthonormalize([ns_vec.vector])

nullspace = PETSc.NullSpace().create(vectors=[ns_vec.vector])

assert nullspace.test(A)

A.setNearNullSpace(nullspace)

Thanks a lot. Why it was not enough to put only V in the locate_dofs_topological function?

Because you are using a mixed-finite element, while the space V is disconnected from the mixed space W.
They relate through the map V_to_W, but this information is not present inside the locate_dofs_topological, thus you end up setting the boundary condition on the wrong dof indices.

1 Like

Sorry to bother, I have another question that is related to the convergence of the problem. This formulation lead me to a convergence error which does not follow any logic: for example I expect that if I halve the mesh also the H¹ norm error should halve for the velocity, but this does not happen. Do you have any suggestion about this aspect?

Please present your code that has the convergence study included.
I’ve shown in other examples how to use ufl to set up the analytical solution. See for instance: Poisson problem with Neumann boundary condition - #6 by dokken
or
Error control: Computing convergence rates — FEniCSx tutorial

Here is the complete code

import numpy as np
import dolfinx
import ufl
from petsc4py import PETSc
from mpi4py import MPI


N = 16
mesh = dolfinx.mesh.create_unit_cube(MPI.COMM_WORLD, N, N, N,
                                     cell_type=dolfinx.mesh.CellType.hexahedron)
P2 = ufl.VectorElement("Lagrange", mesh.ufl_cell(), 2)
P1 = ufl.FiniteElement("Lagrange", mesh.ufl_cell(), 1)
TH = ufl.MixedElement([P2, P1])
W = dolfinx.fem.FunctionSpace(mesh, TH)

def u_ex(x): 
    return np.vstack((4*(x[1]-1)*x[1]*(x[0]-1)*(x[0]-1)*x[0]*x[0]*(x[1]-0.5), -4*(x[0]-0.5)*x[1]*x[1]*x[0]*(x[1]-1)*(x[1]-1)*(x[0]-1), 0*x[0]))
def p_ex(x):
    return x[0]+x[1]+x[2]-1.5

V, V_to_W = W.sub(0).collapse()
Q, Q_to_W = W.sub(1).collapse()

########### NOTE ###############
# Because you are using a mixed-finite element, while the space V is disconnected from the mixed space W.
# They relate through the map V_to_W, but this information is not present inside the locate_dofs_topological, 
# thus you end up setting the boundary condition on the wrong dof indices.

############################### 



### Define the boundary condition -> which in this case are only of Dirichlet type 
uD = dolfinx.fem.Function(V)
uD.interpolate(u_ex)
# Create facet to cell connectivity required to determine boundary facets
tdim = mesh.topology.dim
fdim = tdim - 1
mesh.topology.create_connectivity(fdim, tdim)
boundary_facets = dolfinx.mesh.exterior_facet_indices(mesh.topology)
boundary_dofs = dolfinx.fem.locate_dofs_topological((W.sub(0),V), fdim, boundary_facets)
bcs = [dolfinx.fem.dirichletbc(uD, boundary_dofs, W.sub(0))]

(u, p) = ufl.TrialFunctions(W)
(v, q) = ufl.TestFunctions(W)

x = ufl.SpatialCoordinate(mesh)
f = ufl.as_vector([1-24*(-1 + x[0])*(x[0]-1)*x[0]*x[0]*(-0.5 + x[1])- 48*(0.166667 - x[0] + x[0]*x[0])*(-1 + x[1])*(-0.5 + x[1])*x[1],
                1+24*(-0.5 + x[0])*(-1 + x[1])*(x[1]-1)*x[1]*x[1] + 8*(-1 + x[0])*(-0.5 + x[0])*x[0]*(1 - 6*x[1] + 6*x[1]*x[1]),
                1.0+0*x[0]])

dx = ufl.Measure("dx", domain=mesh)
a = ufl.inner(ufl.grad(u), ufl.grad(v)) * dx - \
    ufl.div(v)*p*dx - q*ufl.div(u)*dx
L = ufl.inner(f, v) * dx

a_f = dolfinx.fem.form(a)
L_f = dolfinx.fem.form(L)
A = dolfinx.fem.petsc.assemble_matrix(a_f, bcs=bcs)
A.assemble()

# Tp print the matrix
'''from scipy.sparse import csr_matrix
assert isinstance(A, PETSc.Mat)
ai, aj, av = A.getValuesCSR()
Asp = csr_matrix((av, aj, ai))'''
#print(Asp)

b = dolfinx.fem.petsc.assemble_vector(L_f)
dolfinx.fem.petsc.apply_lifting(b, [a_f], [bcs])
b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
dolfinx.fem.petsc.set_bc(b, bcs)
b.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)

ns_vec = dolfinx.fem.Function(W)
ns_vec.x.array[Q_to_W] = 1
dolfinx.la.orthonormalize([ns_vec.vector])

nullspace = PETSc.NullSpace().create(vectors=[ns_vec.vector])

assert nullspace.test(A)

A.setNearNullSpace(nullspace)

# Preconditioner matrix form
Pf = dolfinx.fem.form(ufl.inner(ufl.grad(u), ufl.grad(v)) * dx + p * q * dx)
P = dolfinx.fem.petsc.assemble_matrix(Pf, bcs=bcs)
P.assemble()

ksp = PETSc.KSP().create(mesh.comm)
ksp.setOperators(A, P)
ksp.setType(PETSc.KSP.Type.MINRES)
pc = ksp.getPC()
pc.setType("hypre")
pc.setHYPREType("boomeramg")
pc.setUp()

# To view Preconditioner info
# pc.view()


wh = dolfinx.fem.Function(W)
ksp.solve(b, wh.vector)

# To view solver info
# ksp.view()

uh = wh.sub(0).collapse()
p = wh.sub(1).collapse()

print(f"Converged Reason {ksp.getConvergedReason()}" +
      f"\nNumber of iterations {ksp.getIterationNumber()}")

from dolfinx.io import XDMFFile
with XDMFFile(mesh.comm, "stokes/stoke.xdmf", "w") as file:
    file.write_mesh(mesh)
    file.write_function(uh)

with XDMFFile(mesh.comm, "stokes/p.xdmf", "w") as file:
        file.write_mesh(mesh)
        file.write_function(p)

L2_error = dolfinx.fem.form(ufl.inner(uh - uD, uh - uD) * ufl.dx)
error_local = dolfinx.fem.assemble_scalar(L2_error)
error_L2_u = np.sqrt(mesh.comm.allreduce(error_local, op=MPI.SUM))

if mesh.comm.rank == 0:
    print(f"Error_L2 : {error_L2_u:.2e}")

eh = uh - uD
error_H1 = dolfinx.fem.form(ufl.inner(ufl.grad(eh), ufl.grad(eh)) * ufl.dx + ufl.inner(uh - uD, uh - uD) * ufl.dx)
E_H1_u = np.sqrt(mesh.comm.allreduce(dolfinx.fem.assemble_scalar(error_H1), op=MPI.SUM))
if mesh.comm.rank == 0:
    print(f"Error_H1 : {E_H1_u:.2e}")

    
pex = dolfinx.fem.Function(Q)
pex.interpolate(p_ex)

L2_error_p = dolfinx.fem.form((p - pex) * (p - pex) * ufl.dx)
error_local_p = dolfinx.fem.assemble_scalar(L2_error_p)
error_L2_p = np.sqrt(mesh.comm.allreduce(error_local_p, op=MPI.SUM))
if mesh.comm.rank == 0:
    print(f"Error_L2_p : {error_L2_p:.2e}")

I’m not sure what you are talking about here, as the error is clearly dependent on the mesh resolution (and it going down with mesh resolution)
ADding these lines to your script:

h_max = max(dolfinx.cpp.mesh.h(mesh, mesh.topology.dim, np.arange(mesh.topology.index_map(mesh.topology.dim).size_local, dtype=np.int32)))
h_max_global = mesh.comm.allreduce(h_max, op=MPI.MAX)
if mesh.comm.rank == 0:
    print(f"hmax: {h_max_global}")

gives
(N=8)

hmax: 0.21650635094610965
Error_L2 : 1.17e-06
Error_H1 : 2.44e-05
Error_L2_p : 1.36e-03

(N=16)

hmax: 0.10825317547305482
Error_L2 : 1.79e-07
Error_H1 : 6.21e-06
Error_L2_p : 5.61e-04

giving rates of:
u_L2: 2.8
u_H1: 1.97
u_L2_p: 1.27

1 Like

Hi Dokken, if you refine more the grid like N=32,64 the error increase, which is unusual. Those are the results for N = 4,8,16,32,64:

Converged Reason 2
Number of iterations 34
Error_L2 : 1.89e-05
Error_H1 : 2.90e-04
Error_L2_p : 3.89e-03
h: 2.50e-01 Error L2 velocity: 1.89e-05
h: 2.50e-01 Error H1 velocity: 2.90e-04
h: 2.50e-01 Error L2 pressure: 3.89e-03
Converged Reason 2
Number of iterations 35
Error_L2 : 1.17e-06
Error_H1 : 2.44e-05
Error_L2_p : 1.36e-03
h: 1.25e-01 Error L2 velocity: 1.17e-06
h: 1.25e-01 Error H1 velocity: 2.44e-05
h: 1.25e-01 Error L2 pressure: 1.36e-03
Converged Reason 2
Number of iterations 35
Error_L2 : 1.79e-07
Error_H1 : 6.21e-06
Error_L2_p : 5.61e-04
h: 6.25e-02 Error L2 velocity: 1.79e-07
h: 6.25e-02 Error H1 velocity: 6.21e-06
h: 6.25e-02 Error L2 pressure: 5.61e-04
Converged Reason 2
Number of iterations 33
Error_L2 : 1.74e-07
Error_H1 : 1.32e-05
Error_L2_p : 2.55e-04
h: 3.12e-02 Error L2 velocity: 1.74e-07
h: 3.12e-02 Error H1 velocity: 1.32e-05
h: 3.12e-02 Error L2 pressure: 2.55e-04
Converged Reason 2
Number of iterations 31
Error_L2 : 2.25e-07
Error_H1 : 2.55e-05
Error_L2_p : 1.42e-04
h: 1.56e-02 Error L2 velocity: 2.25e-07
h: 1.56e-02 Error H1 velocity: 2.55e-05
h: 1.56e-02 Error L2 pressure: 1.42e-04
Rates L2 velocity: [ 4.00482609  2.71297159  0.04236473 -0.3689221 ]
Rates H1 velocity: [ 3.56969089  1.97319413 -1.08786769 -0.9491815 ]
Rates L2 pressure: [1.51609805 1.27829202 1.13528079 0.84256529]

Since you use an iterative solver, you need to adjust the atol and rtol to make sure you solve the problem sufficiently accurately.

I would start by using a direct solver to eliminate all possible sources of errors (as the solution is then accurate to machine tolerance).

1 Like

Just for curiosity, if instead of the null space I fix the dofs for the pressure in one point to resolve the ill-posedness of the problem, are the following commands right?

p_dof = dolfinx.fem.Function(Q)
with p_dof.vector.localForm() as p_dof_local:
     p_dof_local.set(-1.5)
dofs = dolfinx.fem.locate_dofs_geometrical((W.sub(1),Q),
                    lambda x: np.isclose(x.T, [0, 0, 0]).all(axis=1))
bc_point = dolfinx.fem.dirichletbc(p_dof, dofs, W.sub(1))

In this case I expect that the pressure is fixed on the edge (0,0,0)

The PETSc setNearNullSpace function is for algebraic multigrid solvers, e.g. smoothed aggregation, that require a near nullspace (see PETSc docs).

To eliminate a nullspace, you need to use setNearNullSpace and then a Krylov subspace solver. The preconditioner needs to be able to handle a singular operator. The RHS should be consistent, i.e. not excite the null space.

For LU factorisation, MUMPS had a feature for handling singular systems. See the line https://github.com/FEniCS/dolfinx/blob/deac4f1cba124f09302f8bb501777f3e25a827c4/python/demo/demo_stokes.py#L427.

3 Likes

Good morning,

following the advice of @dokken I saw what happens in the case of a direct solver and everything seams nice until the mesh is not too fine. When a mesh with 64 elements for each edge is set, the algorithm is killing itself without any error message . The code is the following:

# Demo solving Stokes problem in 3D in DOLFINx using MINRES and AMG
# Copyright: 2022
# Author: Jørgen S. Dokken
# License: MIT

import numpy as np
import dolfinx
import ufl
from petsc4py import PETSc
from mpi4py import MPI


def solve_stokes(N):
    mesh = dolfinx.mesh.create_unit_cube(MPI.COMM_WORLD, N, N, N,
                                         cell_type=dolfinx.mesh.CellType.hexahedron)
    P2 = ufl.VectorElement("Lagrange", mesh.ufl_cell(), 2)
    P1 = ufl.FiniteElement("Lagrange", mesh.ufl_cell(), 1)
    TH = ufl.MixedElement([P2, P1])
    W = dolfinx.fem.FunctionSpace(mesh, TH)

    def u_ex(x): 
        return np.vstack((4*(x[1]-1)*x[1]*(x[0]-1)*(x[0]-1)*x[0]*x[0]*(x[1]-0.5), -4*(x[0]-0.5)*x[1]*x[1]*x[0]*(x[1]-1)*(x[1]-1)*(x[0]-1), 0*x[0]))
    def p_ex(x):
        return x[0]+x[1]+x[2]-1.5

    def force(x):
        return np.vstack((1-24*(-1 + x[0])*(x[0]-1)*x[0]*x[0]*(-0.5 + x[1])- 48*(0.166667 - x[0] + x[0]*x[0])*(-1 + x[1])*(-0.5 + x[1])*x[1],
            1+24*(-0.5 + x[0])*(-1 + x[1])*(x[1]-1)*x[1]*x[1] + 8*(-1 + x[0])*(-0.5 + x[0])*x[0]*(1 - 6*x[1] + 6*x[1]*x[1]),
            np.ones_like(x[0])))

    V, V_to_W = W.sub(0).collapse()
    Q, Q_to_W = W.sub(1).collapse()

    ########### NOTE ###############
    # Because you are using a mixed-finite element, while the space V is disconnected from the mixed space W.
    # They relate through the map V_to_W, but this information is not present inside the locate_dofs_topological, 
    # thus you end up setting the boundary condition on the wrong dof indices.

    ############################### 



    ### Define the boundary condition -> which in this case are only of Dirichlet type 
    uD = dolfinx.fem.Function(V)
    uD.interpolate(u_ex)
    # Create facet to cell connectivity required to determine boundary facets
    tdim = mesh.topology.dim
    fdim = tdim - 1
    mesh.topology.create_connectivity(fdim, tdim)
    boundary_facets = dolfinx.mesh.exterior_facet_indices(mesh.topology)
    boundary_dofs = dolfinx.fem.locate_dofs_topological((W.sub(0),V), fdim, boundary_facets)
    bc1 = dolfinx.fem.dirichletbc(uD, boundary_dofs, W.sub(0))

    zero = dolfinx.fem.Function(Q)
    zero.x.set(-1.5)
    dofs = dolfinx.fem.locate_dofs_geometrical((W.sub(1), Q), lambda x: np.isclose(x.T, [0, 0, 0]).all(axis=1))
    bc2 = dolfinx.fem.dirichletbc(zero, dofs, W.sub(1))

    bcs = [bc1, bc2]

    (u, p) = ufl.TrialFunctions(W)
    (v, q) = ufl.TestFunctions(W)

    x = ufl.SpatialCoordinate(mesh)
    f = ufl.as_vector([1-24*(-1 + x[0])*(x[0]-1)*x[0]*x[0]*(-0.5 + x[1])- 48*(0.166667 - x[0] + x[0]*x[0])*(-1 + x[1])*(-0.5 + x[1])*x[1],
                    1+24*(-0.5 + x[0])*(-1 + x[1])*(x[1]-1)*x[1]*x[1] + 8*(-1 + x[0])*(-0.5 + x[0])*x[0]*(1 - 6*x[1] + 6*x[1]*x[1]),
                    1.0+0*x[0]])
    ff = dolfinx.fem.Function(V)
    ff.interpolate(force)
    dx = ufl.Measure("dx", domain=mesh)
    a = ufl.inner(ufl.grad(u), ufl.grad(v)) * dx - \
        ufl.div(v)*p*dx - q*ufl.div(u)*dx
    L = ufl.inner(ff, v) * dx

    a_f = dolfinx.fem.form(a)
    L_f = dolfinx.fem.form(L)
    A = dolfinx.fem.petsc.assemble_matrix(a_f, bcs=bcs)
    A.assemble()

    # Tp print the matrix
    '''from scipy.sparse import csr_matrix
    assert isinstance(A, PETSc.Mat)
    ai, aj, av = A.getValuesCSR()
    Asp = csr_matrix((av, aj, ai))'''
    #print(Asp)

    b = dolfinx.fem.petsc.assemble_vector(L_f)
    dolfinx.fem.petsc.apply_lifting(b, [a_f], [bcs])
    b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
    dolfinx.fem.petsc.set_bc(b, bcs)
    #b.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)


    # Create and configure solver
    ksp = PETSc.KSP().create(mesh.comm)
    ksp.setOperators(A)
    ksp.setType("preonly")
    ksp.getPC().setType("lu")
    #ksp.getPC().setFactorSolverType("superlu_dist")
    ksp.getPC().setFactorSolverType("mumps")
    opts = PETSc.Options()
    opts["mat_mumps_icntl_14"] = 80  # Increase MUMPS working memory
    opts["mat_mumps_icntl_24"] = 0  # Option to support solving a singular matrix (pressure nullspace)
    opts["mat_mumps_icntl_25"] = 0  # Option to support solving a singular matrix (pressure nullspace)
    opts["ksp_error_if_not_converged"] = 1
    ksp.setFromOptions()

    # Compute the solution
    U = dolfinx.fem.Function(W)
    ksp.solve(b, U.vector)

    # Split the mixed solution and collapse
    uh = U.sub(0).collapse()
    p = U.sub(1).collapse()

    print(f"Converged Reason {ksp.getConvergedReason()}" +
          f"\nNumber of iterations {ksp.getIterationNumber()}")

    from dolfinx.io import XDMFFile
    with XDMFFile(mesh.comm, "stokes/stoke.xdmf", "w") as file:
        file.write_mesh(mesh)
        file.write_function(uh)

    with XDMFFile(mesh.comm, "stokes/p.xdmf", "w") as file:
            file.write_mesh(mesh)
            file.write_function(p)

    L2_error = dolfinx.fem.form(ufl.inner(uh - uD, uh - uD) * ufl.dx)
    error_local = dolfinx.fem.assemble_scalar(L2_error)
    error_L2_u = np.sqrt(mesh.comm.allreduce(error_local, op=MPI.SUM))

    if mesh.comm.rank == 0:
        print(f"Error_L2 : {error_L2_u:.2e}")

    eh = uh - uD
    error_H1 = dolfinx.fem.form(ufl.inner(ufl.grad(eh), ufl.grad(eh)) * ufl.dx + ufl.inner(uh - uD, uh - uD) * ufl.dx)
    E_H1_u = np.sqrt(mesh.comm.allreduce(dolfinx.fem.assemble_scalar(error_H1), op=MPI.SUM))
    if mesh.comm.rank == 0:
        print(f"Error_H1 : {E_H1_u:.2e}")

        
    pex = dolfinx.fem.Function(Q)
    pex.interpolate(p_ex)

    L2_error_p = dolfinx.fem.form((p - pex) * (p - pex) * ufl.dx)
    error_local_p = dolfinx.fem.assemble_scalar(L2_error_p)
    error_L2_p = np.sqrt(mesh.comm.allreduce(error_local_p, op=MPI.SUM))
    if mesh.comm.rank == 0:
        print(f"Error_L2_p : {error_L2_p:.2e}")

    return error_L2_u, E_H1_u, error_L2_p



Ns = [8, 16, 32, 64]
L2_u = np.zeros(len(Ns), dtype=PETSc.ScalarType)
H1_u = np.zeros(len(Ns), dtype=PETSc.ScalarType)
L2_p = np.zeros(len(Ns), dtype=PETSc.ScalarType)
hs = np.zeros(len(Ns), dtype=np.float64)
for i, N in enumerate(Ns):
    error_L2_u, E_H1_u, error_L2_p = solve_stokes(N)

    L2_u[i] = error_L2_u
    H1_u[i] = E_H1_u
    L2_p[i] = error_L2_p
    hs[i] = 1./Ns[i]
    if MPI.COMM_WORLD.rank == 0:
        print(f"h: {hs[i]:.2e} Error L2 velocity: {L2_u[i]:.2e}")
        print(f"h: {hs[i]:.2e} Error H1 velocity: {H1_u[i]:.2e}")
        print(f"h: {hs[i]:.2e} Error L2 pressure: {L2_p[i]:.2e}")

rates_L2u = np.log(L2_u[1:]/L2_u[:-1])/np.log(hs[1:]/hs[:-1])
rates_H1u = np.log(H1_u[1:]/H1_u[:-1])/np.log(hs[1:]/hs[:-1])
rates_L2p = np.log(L2_p[1:]/L2_p[:-1])/np.log(hs[1:]/hs[:-1])
if MPI.COMM_WORLD.rank == 0:
    print(f"Rates L2 velocity: {rates_L2u}")
    print(f"Rates H1 velocity: {rates_H1u}")
    print(f"Rates L2 pressure: {rates_L2p}")

Concerning the output:

Converged Reason 4
Number of iterations 1
Error_L2 : 1.18e-06
Error_H1 : 2.21e-05
Error_L2_p : 2.70e-05
h: 1.25e-01 Error L2 velocity: 1.18e-06
h: 1.25e-01 Error H1 velocity: 2.21e-05
h: 1.25e-01 Error L2 pressure: 2.70e-05
Converged Reason 4
Number of iterations 1
Error_L2 : 7.42e-08
Error_H1 : 2.06e-06
Error_L2_p : 1.97e-06
h: 6.25e-02 Error L2 velocity: 7.42e-08
h: 6.25e-02 Error H1 velocity: 2.06e-06
h: 6.25e-02 Error L2 pressure: 1.97e-06
Converged Reason 4
Number of iterations 1
Error_L2 : 7.41e-09
Error_H1 : 2.16e-07
Error_L2_p : 2.51e-07
h: 3.12e-02 Error L2 velocity: 7.41e-09
h: 3.12e-02 Error H1 velocity: 2.16e-07
h: 3.12e-02 Error L2 pressure: 2.51e-07
Killed

Do you have any suggestion about what can I do wrong? Thanks in advance

That usually means that you are running out of memory on your system

The strange fact is that if I run the same algorithm with iterative solver, but to remove the singularity I fix the value of the pressure in one node, instead of removing the null space. Then the algorithm gave results which are not good (meaning that the error for the pressure tend to increase for the case of 64 elements for each edge), but it is able to run for the case of a more fine grid.