# Pressure error is too large

**URL:** <https://fenicsproject.discourse.group/t/pressure-error-is-too-large/3610>\
**Category:** variational formulation\
**Created:** [June 22, 2020, 5:32pm UTC](https://fenicsproject.discourse.group/t/pressure-error-is-too-large/3610 "2020-06-22T17:32:39Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![Khan](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/khan/32/1399_2.png) [@Khan](https://fenicsproject.discourse.group/u/Khan)\
**Post date:** [June 22, 2020, 5:32pm UTC](https://fenicsproject.discourse.group/t/pressure-error-is-too-large/3610/1 "2020-06-22T17:32:39Z")

</div>

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(pi_x[0])  
value[1] = -pi_x[1]_cos(pi_x[0])  
value[2] = sin(pi\*x[0])_cos(pi_x[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(pi_x[0])  
value[1] = -pi_x[1]_cos(pi_x[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(pi_x[0]);  
u1x = pi_cos(pi_x[0]);  
u1y = 0;  
u1lap = -pi_pi_sin(pi_x[0]);  
u2 = -pi_x[1]cos(pix[0]);  
u2x = pi_pi_x[1]sin(pix[0]);  
u2y = -pi_cos(pi_x[0]);  
u2lap = pi_pi_pi_x[1]_cos(pi_x[0]);  
px = pi_cos(pi_x[0])_cos(pi_x[1]);  
py = -pi_sin(pi_x[0])_sin(pi_x[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(pi_x[0])  
value[1] = -pi_x[1]_cos(pi_x[0])  
value[2] = sin(pi\*x[0])_cos(pi_x[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)\*\*2_dx  
ep = (p-p\_e)\*\*2_dx  
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 + q_div(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

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [June 23, 2020, 12:16pm UTC](https://fenicsproject.discourse.group/t/pressure-error-is-too-large/3610/2 "2020-06-23T12:16:30Z")

</div>

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:

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

```

---

<div class="post-metadata">

**Author:** ![Khan](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/khan/32/1399_2.png) [@Khan](https://fenicsproject.discourse.group/u/Khan)\
**Post date:** [June 23, 2020, 1:48pm UTC](https://fenicsproject.discourse.group/t/pressure-error-is-too-large/3610/3 "2020-06-23T13:48:42Z")

</div>

Thanks for your explanation. It works.  
Best regards  
Khan