Convergence in Stokes problem

Good evening,

in the Stokes problem using P1-P1 finite elements with a stabilization term (I implemented the easiest, i.e. the Brezzi-Pitkaranta stabilization), one should expect a convergence like O(h) for the H1 error in the velocity and the pressure. In my case this does not happens and following previous advice I tried first to solve this problem using a direct solver, but still the results are not right.
The minimal code is the following:

import numpy as np

import ufl
from dolfinx import cpp as _cpp
from dolfinx import fem
from dolfinx.fem import (Constant, Function, FunctionSpace, dirichletbc,
                         extract_function_spaces, form,
                         locate_dofs_geometrical, locate_dofs_topological)
from dolfinx.io import XDMFFile
from dolfinx.mesh import (CellType, GhostMode, create_rectangle, create_unit_cube,
                          locate_entities_boundary, exterior_facet_indices)
from ufl import div, dx, grad, inner

from mpi4py import MPI
from petsc4py import PETSc


# Create mesh
N = 32
msh = create_unit_cube(MPI.COMM_WORLD,
                       N, N, N,
                       CellType.hexahedron, GhostMode.none)


'''# Function to mark x = 0, x = 1 and y = 0
def noslip_boundary(x):
    return np.logical_or(np.logical_or(np.isclose(x[0], 0.0),
                                       np.isclose(x[0], 1.0)),
                         np.isclose(x[1], 0.0))


# Function to mark the lid (y = 1)
def lid(x):
    return np.isclose(x[1], 1.0)


# Lid velocity
def lid_velocity_expression(x):
    return np.stack((np.ones(x.shape[1]), np.zeros(x.shape[1])))'''

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]-1.5

def f_ex(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])))


P2 = ufl.VectorElement("Lagrange", msh.ufl_cell(), 1)
P1 = ufl.FiniteElement("Lagrange", msh.ufl_cell(), 1)

# Create the function space
TH = P2 * P1
W = FunctionSpace(msh, TH)
W0, _ = W.sub(0).collapse()
W1, _ = W.sub(1).collapse()

'''# No slip boundary condition
noslip = Function(V)
facets = locate_entities_boundary(msh, 1, noslip_boundary)
dofs = locate_dofs_topological((W.sub(0), V), 1, facets)
bc0 = dirichletbc(noslip, dofs, W.sub(0))


# Driving velocity condition u = (1, 0) on top boundary (y = 1)
lid_velocity = Function(W0)
lid_velocity.interpolate(lid_velocity_expression)
facets = locate_entities_boundary(msh, 1, lid)
dofs = locate_dofs_topological((W.sub(0), V), 1, facets)
bc1 = dirichletbc(lid_velocity, dofs, W.sub(0))'''

uD = Function(W0)
uD.interpolate(u_ex)
uD.vector.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
# Create facet to cell connectivity required to determine boundary facets
tdim = msh.topology.dim
fdim = tdim - 1
msh.topology.create_connectivity(fdim, tdim)
boundary_facets = exterior_facet_indices(msh.topology)
boundary_dofs = locate_dofs_topological((W.sub(0),W0), fdim, boundary_facets)
bc1 = dirichletbc(uD, boundary_dofs, W.sub(0))

# Since for this problem the pressure is only determined up to a
# constant, we pin the pressure at the point (0, 0)
zero = Function(W1)
zero.x.set(-1.5)
zero.vector.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
dofs = locate_dofs_geometrical((W.sub(1), W1), lambda x: np.isclose(x.T, [0, 0, 0]).all(axis=1))
bc2 = dirichletbc(zero, dofs, W.sub(1))

# Collect Dirichlet boundary conditions
bcs = [bc1, bc2]

# Define variational problem
(u, p) = ufl.TrialFunctions(W)
(v, q) = ufl.TestFunctions(W)
f = Function(W0)
f.interpolate(f_ex)
h = PETSc.ScalarType(ufl.sqrt(2)*1./N)
beta = PETSc.ScalarType(0.1)
a = form((inner(grad(u), grad(v)) + inner(p, div(v)) + inner(div(u), q)) * dx - h*h*beta*ufl.inner(ufl.nabla_grad(p),ufl.nabla_grad(q))*dx)
L = form(inner(f, v) * dx)


# Assemble LHS matrix and RHS vector
A = fem.petsc.assemble_matrix(a, bcs=bcs)
A.assemble()
b = fem.petsc.assemble_vector(L)

fem.petsc.apply_lifting(b, [a], bcs=[bcs])
b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)

# Set Dirichlet boundary condition values in the RHS
fem.petsc.set_bc(b, bcs)

# Create and configure solver
ksp = PETSc.KSP().create(msh.comm)
ksp.setOperators(A)
ksp.setType("preonly")
ksp.getPC().setType("lu")
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
opts["ksp_monitor"] = None
ksp.setFromOptions()

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

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

L2_error = fem.form(ufl.inner(u - uD, u - uD) * ufl.dx)
error_local = fem.assemble_scalar(L2_error)
error_L2 = np.sqrt(msh.comm.allreduce(error_local, op=MPI.SUM))

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

eh = u - uD
error_H10 = fem.form(ufl.inner(ufl.grad(eh), ufl.grad(eh)) * ufl.dx)
E_H10 = np.sqrt(msh.comm.allreduce(fem.assemble_scalar(error_H10), op=MPI.SUM))
if msh.comm.rank == 0:
    print(f"H01-error: {E_H10:.2e}")

pD = Function(W1)
pD.interpolate(p_ex)
L2_errorp = fem.form(ufl.inner(p - pD, p - pD) * ufl.dx)
error_localp = fem.assemble_scalar(L2_errorp)
error_L2p = np.sqrt(msh.comm.allreduce(error_localp, op=MPI.SUM))

if msh.comm.rank == 0:
    print(f"Error_L2p : {error_L2p:.2e}")

# Write the solution to file
with XDMFFile(MPI.COMM_WORLD, "out_stokes/u.xdmf", "w") as ufile_xdmf:
    ufile_xdmf.write_mesh(msh)
    ufile_xdmf.write_function(u)

with XDMFFile(MPI.COMM_WORLD, "out_stokes/p.xdmf", "w") as pfile_xdmf:
    pfile_xdmf.write_mesh(msh)
    pfile_xdmf.write_function(p)

with XDMFFile(MPI.COMM_WORLD, "out_stokes/uex.xdmf", "w") as pfile_xdmf:
    pfile_xdmf.write_mesh(msh)
    pfile_xdmf.write_function(uD)

The error for a mesh with respectively 8, 16, 32 elements are:

Error_L2 : 5.08e-03
H01-error: 3.39e-02
Error_L2p : 2.67e+00

Error_L2 : 1.41e-03
H01-error: 1.14e-02
Error_L2p : 2.90e+00

Error_L2 : 3.68e-04
H01-error: 3.71e-03
Error_L2p : 3.03e+00

Which are note corrects (in particular the one of the pressure??). Can someone give me an hand on this problem that is driving me crazy?

Looks like you have the wrong sign in front of the pressure term.

a = form((inner(grad(u), grad(v)) + inner(p, div(v)) + inner(div(u), q)) * dx - h*h*beta*ufl.inner(ufl.nabla_grad(p),ufl.nabla_grad(q))*dx)

You could either flip the signs on the continuity and pressure terms to maintain symmetry, or just post process your pressure solution to account for it.

Thanks for the answer. This solved the bed convergence of the pressure, but still I don’t get the rate of convergence I expect.
Did you notice something else that is wrong?