Hello,
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)
hvalues.append(mesh.hmin())
# 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,
with
\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.
