Stokes problem (pressure order convergence)

Dear developer,
I am new to NGSolve and trying to solve the Stokes problem to check the velocity and pressure convergence order for the known solution. The velocity error converges to the right order but the pressure doesn't converge. Anyone can help me to get the right convergence orders please? Thanks
Best regards,
Martin

The code:

from ngsolve import *
# import netgen.gui
from netgen.geom2d import unit_square
from ngsolve.webgui import Draw
from math import pi
from netgen.geom2d import SplineGeometry
import numpy as np

geo = SplineGeometry()
geo.AddRectangle((0,0),(1,1),bc="rectangle")

nu = 1.0
# exact velocity
ex_vel = (sin(2*pi*x) * sin(2*pi*y),
          cos(2*pi*x) * cos(2*pi*y))
# exact pressure
ex_pre = 0.25*(cos(4*pi*x) - cos(4*pi*y))
# source term
src = ( 8*pi*pi*sin(2*pi*x)*sin(2*pi*y) + pi*sin(4*pi*x),
       8*pi*pi*cos(2*pi*y)*cos(2*pi*x) + pi*sin(4*pi*y))

mesh = Mesh(geo.GenerateMesh(maxh=1/(2*N)))
mesh.Curve(order=order)
print("ElementBoundary=",mesh.GetBoundaries())
print("NumEles=",mesh.ne)
print("NumEdges=",mesh.nedge)
print("NumVertex=",mesh.nv)
    
# finite element formulation
V = VectorH1(mesh,order=order, dirichlet="rectangle")
Q = H1(mesh,order=order-1)
    
X = FESpace([V,Q])

u,p = X.TrialFunction()
v,q = X.TestFunction()

print("velocity dofs ", V.ndof, X.Range(0))
print("pressure dofs ", Q.ndof, X.Range(1))
print("Total dofs ", X.ndof)
    
# assembling the system matrix
stokes = nu*InnerProduct(grad(u), grad(v))-div(u)*q-div(v)*p - 1e-10*p*q
a = BilinearForm(X)
a += SymbolicBFI(stokes)
a.Assemble()
# assembling the rhs
f = LinearForm(X)
source = CoefficientFunction(src)
f += SymbolicLFI(InnerProduct(source,v))
f.Assemble()
    
#Exact Dirichelt BC for velocity
gfu = GridFunction(X)

u_exact = CoefficientFunction(ex_vel)
gfu.components[0].Set(u_exact, definedon=mesh.Boundaries("rectangle"))    
    
# solve the linear system
inv_stokes = a.mat.Inverse(X.FreeDofs())

res = f.vec.CreateVector()
res.data = f.vec - a.mat*gfu.vec
gfu.vec.data += inv_stokes * res

velocity = gfu.components[0]
pressure = gfu.components[1]
    
# Draw( velocity, mesh, "velocity" )
# Draw( Norm(velocity), mesh, "absvel(hdiv)")
# Draw( pressure, mesh, "pressure" )
    
# computing the error estimates
p_exact = CoefficientFunction(ex_pre)
pressure = CoefficientFunction(gfu.components[1])
errp = sqrt(Integrate((p_exact-pressure)*(p_exact-pressure), mesh))

velocity = CoefficientFunction(gfu.components[0])    
erru = sqrt (Integrate ( (u_exact-velocity )*(u_exact-velocity), mesh))

print("L2(u):", erru)
print("L2(p):", errp) 

Hi,

your source term has a small typo. The sign in the pressure part should be negative, i.e.

src = ( 8*pi*pi*sin(2*pi*x)*sin(2*pi*y) - pi*sin(4*pi*x),
       8*pi*pi*cos(2*pi*y)*cos(2*pi*x) + pi*sin(4*pi*y))

Best
Tim

Dear Tim,
thank you very much. It was a stupid mistake but I appreciate your help.
Best regards,
Martin