Why does pressure error do not decrease as mesh is refined? (Stokes flow)

I am trying to generate convergence plots for my project, and I have noticed that my pressure error does not decrease with decreasing mesh size, whereas the velocity behaves as expected. Does this mean that I am not converging to the right pressure value? I know that depending on the combination of boundary conditions used, pressure is unique up to a constant. My domain is a rectangle, but the top is split into two halves. Only the top right boundary has a Neumann zero stress condition, which I would have thought would fix the pressure to a unique value.

This is a minimum working example where I print the resulting error values:

from dolfin import *
import numpy as np

def epsilon(v):
    return sym(as_tensor([[v[0].dx(0), v[0].dx(1), 0],
                          [v[1].dx(0), v[1].dx(1), 0],
                          [0, 0, 0]]))

#stress tensor
def sigma(v, p):
    return 2*mu*epsilon(v)-p*Identity(3)

def cond(v):
    return sym(as_tensor([[v[0].dx(0), v[0].dx(1)],
                          [v[1].dx(0), v[1].dx(1)]]))

def sigmabc(v, p):
    return 2*mu*cond(v) - p*Identity(2)

#Initial coarse mesh
mesh = RectangleMesh(Point(0, 0), Point(0.2, 1.0), 3, 3)

#Manufactured solution of the easy case
solns_ex = ['0.0', '-(1-x[1])*(1-x[1])*(x[0]*x[0]-0.04)', '0.0', '2*(1-x[1])*(x[0]*x[0]-0.04)']
solns_f = ['0.0', '2*(x[0]*x[0]-0.04) + 2*(1-x[1])*(1-x[1])']

mu = Constant(1.0) # constant viscosity

inflow = 'near(x[1], 1.0) && x[0]<=0.1'
right = 'near(x[1], 1.0) && x[0]>=0.1'
wall = 'near(x[0], 0.2)'
centre = 'near(x[0], 0.0)'
outflow = 'near(x[1], 0.0)'

#initialise variables
hvalues = []
errors_u = []
errors_p = []
iterations = 6 # total number of iterations

for i in range(iterations):

    n = FacetNormal(mesh)

    # Define Taylor--Hood function space W
    V = VectorElement("CG", triangle, 2)  # original spaces
    Q = FiniteElement("CG", triangle, 1)  # original spaces
    W = FunctionSpace(mesh, MixedElement([V, Q]))

    # Define Function and TestFunction(s)
    w = Function(W)
    (u, p) = split(w)
    (v, q) = split(TestFunction(W))

    # Manufactured solution of the easy case
    u_ex = Expression((solns_ex[0], solns_ex[1]), element=V, domain=mesh)
    p_ex = Expression(solns_ex[2], element=Q, domain=mesh)
    g_ex = Expression(solns_ex[3], element=Q, domain=mesh)
    f_ex = Expression((solns_f[0], solns_f[1]), element=V, domain=mesh)

    bcu_inflow = DirichletBC(W.sub(0), u_ex, inflow)
    bcu_wall = DirichletBC(W.sub(0), u_ex, wall)
    bcu_outflow = DirichletBC(W.sub(0), u_ex, outflow)
    bcu_centre = DirichletBC(W.sub(0).sub(0), Constant(0.0), centre) # this is a SYMMETRY condition, we also need to set u[1].dx(0) == 0
    bcs = [bcu_inflow, bcu_wall, bcu_outflow, bcu_centre]

    x = SpatialCoordinate(mesh)

    facet_f = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)  # FACET function
    CompiledSubDomain('near(x[1], 1.0) && x[0]<=0.1').mark(facet_f, 0)
    CompiledSubDomain('near(x[1], 1.0) && x[0]>=0.1').mark(facet_f, 1)
    CompiledSubDomain('near(x[0], 0.2)').mark(facet_f, 2)  # wall
    CompiledSubDomain('near(x[1], 0.0)').mark(facet_f, 3)  # outflow
    CompiledSubDomain('near(x[0], 0.0)').mark(facet_f, 4)

    # Create the measure
    ds = Measure("ds", domain=mesh, subdomain_data=facet_f)
    a1 = (inner(sigma(u, p), epsilon(v))) * dx
    a2 = (- div(u) * q + g_ex * q - dot(f_ex, v)) * dx
    b1 = - dot(dot(sigmabc(u, p), n), v) * ds(0) - dot(dot(sigmabc(u, p), n), v) * ds(2) - dot(dot(sigmabc(u, p), n),
                                                                                               v) * ds(3) \
         + dot(p*n, v) * ds(4) # we include this as the normal component of the viscous stress tensor is zero at ds(4)
    F = a1 + a2 + b1
    # Solve problem
    solve(F == 0, w, bcs)
    (u, p) = w.split()
    errors_u.append(errornorm(u_ex, u))
    errors_p.append(errornorm(p_ex, p))
    mesh = refine(mesh)

print('The u errors are', errors_u)
print('The p errors are', errors_p)

# The stress at ds(1) should be zero.
Vsig = TensorFunctionSpace(mesh, "DG", degree=0)
sig = Function(Vsig, name="Stress")
sig.assign(project(sigmabc(u, p), Vsig))
area1 = assemble(1.0 * ds(1))
normal_stress1 = assemble(inner(sig * n, n) * ds(1))/area1
print("Normal stress at boundary 1 (should be zero)", normal_stress1)

The resulting u errors are [6.80482422516266e-05, 7.112756256197105e-05, 1.754189837557531e-05, 4.3576380166020975e-06, 1.0870749955242577e-06, 2.7153954390694644e-07]
and the p errors are [0.01216340102162744, 0.009082409735152962, 0.012222625830076746, 0.014378000041386134, 0.014911394032883796, 0.015041986916233542]

I think I should expect that p should have an error of \mathcal{O}(10^{-2}) as I am using P1 spaces for it. However, it sometimes increases with mesh refinement which worries me a bit. I don’t know if I can attribute it to standard numerical error.

For reference, the manufactured solution I am using is (the simplest one I could think of):
\mathbf{u} = \left(u, v \right) =\left( 0, -(1-y^2)(x^2-0.04)\right), and p\equiv 0, which satisfy the equations:
\frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} = g, \\ \frac{\partial p}{\partial x} = \frac{\partial ^2 u}{\partial x ^2} + \frac{\partial ^2 u}{\partial y ^ 2} + f_u, \\ \frac{\partial p}{\partial y} = \frac{\partial ^2 v}{\partial x ^2} + \frac{\partial ^2 v}{\partial y ^ 2} + f_v,
\mathbf{f} = \left(f_u, f_v \right) = \left(0, 2(1-y)^2 (x^2 - 0.04) \right), \\ g = 2(1-y)(x^2 -0.04)
and all boundary conditions given by evaluating the exact solution, apart from the symmetry condition, which I impose as
u = 0, \quad \frac{\partial v}{\partial x} = 0 \quad on \quad \partial \Omega_4 (ds(4) in the code)

I have tried with three other manufactured solutions and I get the same behaviour, so it must be something I am misunderstanding about the implementation itself.

Let me know if anything comes to mind, I would really appreciate it.

Some things you could check:

  1. Are the boundary marker facet_f as expected? (Save it to pvd and visualize in Paraview).
  2. Other things to check is how does the solution of the pressure compare to the exact solution visually? (Either by projecting the error into an appropriate function space, or using the “eyeball”-norm.
  3. Out of curiosity, why epsilon and sigma defined as 3x3 matrices, while your mesh is 2D (and would suggest that you should use 2x2 matrices), like epsilon = lambda u: ufl.sym(ufl.grad(u))

If neither of these suggestions gets you any further, I would suggest starting with a simpler problem, have you considered using a manufactured solution for the “incompressible stokes equations” first?

Thank you so much for taking the time to respond to my post.

  1. The facets are as expected.
  2. The pressure solution has its largest error near the bottom of the domain, almost like the outflow condition for u makes the pressure struggle. This is what the error looks like at the most refined level:
  3. I personally struggle (a lot) with changing between coordinate systems. Since I have to run various more complicated examples in cylindrical coordinates and compare to Cartesian, I wanted to prevent future mistakes by always writing 3x3 tensors, as the cylindrical tensor has an extra term. I could also change the variational formulation to include that term, so it’s really the same. I should have made them 2D for the MWE, sorry about that.

I finally ran an incompressible problem, and you’re right. The error for the velocity and the pressure now behaves as needed. I based the manufactured solution on this post in steady state. I was aware that there are some issues for compressible flow, I was focused on other aspects of my solution and didn’t connect the dots.

Thank you!