Pressure error is too large

variational formulation
1

Dear all,
I am new to fenics and tried to solve a simple stokes problem using P2/P1 elements. I have pure Dirichlet boundary conditions and my exact solution looks like
class ExactUserExpression):
def eval(self, value, x):
value[0] = sin(pix[0])
value[1] = -pix[1]cos(pix[0])
value[2] = sin(pi*x[0])cos(pix[1])
def value_shape(self):
return(3,)

I computed the errors after solving the Stokes problem but I am getting a very large pressure error. It seems that I have a problem of integral mean. Can any one guide me how to fix integral mean constraint via Lagrange multiplier or any other way to solve that issue?

from dolfin import *
import matplotlib.pyplot as plt
import numpy as np

Load mesh and subdomains

nx = ny = 4
mesh = UnitSquareMesh(nx, ny)
#plot(mesh, ‘grid 1’)
#plt.show()

P2 = VectorElement(“CG”, mesh.ufl_cell(), 2)
P1 = FiniteElement(“CG”, mesh.ufl_cell(), 1)
TH = P2 * P1
W = FunctionSpace(mesh, TH)

##Dirichlet bc everywhere
def dir_bound(x, on_boundary):
return on_boundary

#bc_p = DirichletBC(W.sub(1), 0., “x[0] < DOLFIN_EPS && x[1] < DOLFIN_EPS”, “pointwise”)

class sc_velocity(UserExpression):
def eval(self, value, x):
value[0] = sin(pix[0])
value[1] = -pix[1]cos(pix[0])
def value_shape(self):
return (2,)

u_D = sc_velocity(degree = 2)
bc_u = DirichletBC(W.sub(0), u_D, dir_bound)
bc = [bc_u]

eps = 1.0
class source_term(UserExpression):
def eval(self, value, x):
u1 = sin(pix[0]);
u1x = picos(pix[0]);
u1y = 0;
u1lap = -pipisin(pix[0]);
u2 = -pix[1]cos(pix[0]);
u2x = pipix[1]sin(pix[0]);
u2y = -picos(pix[0]);
u2lap = pipipix[1]cos(pix[0]);
px = picos(pix[0])cos(pix[1]);
py = -pisin(pix[0])sin(pix[1]);

value[0] = -nu*u1lap + px
value[1] = -nu*u2lap + py

def value_shape(self):
return (2,)

class exact(UserExpression):
def eval(self, value, x):
value[0] = sin(pix[0])
value[1] = -pix[1]cos(pix[0])
value[2] = sin(pi*x[0])cos(pix[1])
def value_shape(self):
return(3,)

b = sc_velocity(degree=2)
f = source_term(degree = 2)

error computation

def compute_errors(up, ex, name):

get function space

W = up.function_space()
V = W.sub(0)
P = W.sub(1)
u, p = up.split()
u_e,p_e = ex.split()
eu = (u-u_e)**2dx
ep = (p-p_e)**2dx
grad_eu = grad(u-u_e)**2*dx
l2u = sqrt(abs(assemble(eu)))
l2p = sqrt(abs(assemble(ep)))
h1 = sqrt(abs(assemble(grad_eu)))
print(name, ‘\nL2(u) = %.6g: H1(u)=%0.6g: L2§=%6g’ % (l2u, h1, l2p))

def test_problem():
print(“Solving plan Scott-Vogelius Galerkin”)
#variational formulation
(u, p) = TrialFunctions(W)
(v, q) = TestFunctions(W)
#bilinear and linear form
a = nu*inner(grad(u), grad(v))dx - div(v)pdx + qdiv(u)*dx
L = inner(f, v)*dx

up = Function(W)
solve(a==L,up,bc)

#error computation
exact = exact(degree = 4)
ex = interpolate(exact, W)
compute_errors(up, ex, “FEM”)
if name== “main”:
nu = Constant(1e-0)
test_problem()

Thanks for your help
Khan

2

For anyone to be able to run your code, you need to encapsulate it with ```

You need to subtract the mean value as the pressure can only be determined up to a constant since you only have Dirichlet conditions on the velocity field. This can be done in the following way:

u, p = up.split(deepcopy=True)
p.vector()[:] -= assemble(p*dx)/assemble(1*dx(domain=mesh))
3

Thanks for your explanation. It works.
Best regards
Khan