How to impose a boundary condition on P0 elements for the Stokes problem?

variational formulation
1

Hello,

I use FEniCS with Docker on Ubuntu 18.04. I am trying to solve the steady Stokes problem on the unit square \Omega = (0,1)² :

-\Delta v+\nabla p = f \quad \text { in } \Omega \\ \nabla \cdot v =0 \quad \text { in } \Omega \\ v =0 \quad \text { on } \partial \Omega

with the prescribed body force f(x,y) = (-4 {\pi}^2 \sin(2 \pi y), -16{\pi}^2 \sin(2 \pi x) \cos(2 \pi y) + 4 {\pi}^2 \sin(2 \pi x)) corresponding to the manufactured solutions

v(x,y) = \begin{pmatrix} (\cos(2 \pi x)-1) \sin(2 \pi y) \\ - \sin(2 \pi x)(\cos(2 \pi y)-1) \end{pmatrix} \text{ and } p(x,y) = -4 \pi \sin(2 \pi x) \sin(2 \pi y).

I want to use P^2 - P^0 mixed finite elements for the velocity v and the pressure p.

Here is my code :

from dolfin import *

# Load mesh 
mesh = UnitSquareMesh(32, 32)

# Define function spaces
Elem_u = VectorElement("CG", mesh.ufl_cell(), 2)
Elem_p = FiniteElement("DG", mesh.ufl_cell(), 0)
W = FunctionSpace(mesh, Elem_u * Elem_p)

# Define Dirichlet boundary
def boundary(x):
    return x[0] < DOLFIN_EPS or x[0] > 1.0 - DOLFIN_EPS or x[1] < DOLFIN_EPS or x[1] > 1.0 - DOLFIN_EPS

# Boundary conditions for velocity and pressure
bc0 = DirichletBC(W.sub(0), Constant((0.0, 0.0)), boundary)        
bc1 = DirichletBC(W.sub(1), Constant(0.0), boundary, "pointwise")
bcs = [bc0, bc1]

# Define variational problem
(u, p) = TrialFunctions(W)
(v, q) = TestFunctions(W)
f = Expression(("-4*pi*pi*sin(2*pi*x[1])","-16*pi*pi*sin(2*pi*x[0])*cos(2*pi*x[1]) + 4*pi*pi*sin(2*pi*x[0])"), degree=2)
a = (inner(grad(u), grad(v)) - div(v)*p + q*div(u))*dx
L = inner(f, v)*dx

# Compute solution
w = Function(W)
solve(a == L, w, bcs)

There is no error message but the pressure obtained does not correspond to the manufactured solution. Indeed, it seems that the zero boundary condition for the pressure (bc1 in the code) is not correctly imposed, as illustrates the following plot :

p

Do you know where the problem could come from ?

Note that if I use P^2 - P^1 elements I get the correct solution, so that I really think the problem comes from imposing a Dirichlet boundary conditions on P^0 elements.

Thank you in advance for your help :slight_smile:

2

An educated guess here is that you are trying to pin the pressure at the boundaries. However, for a DG 0 function, the degree of freedom is at the cell center, Which will never be close to 0 or 1 with Machine tolerance. You Need to create separate boundary markers for the pressure and velocity, were the pressure is checking within a bigger distance than DOLFIN_EPS

2 Likes
3

Shouldn’t the pressure be imposed either by adding the constraint \int_\Omega p \; \mathrm{d}x = 0, or by weak imposition of the BCs? Prescribing BC values to DoFs seems risky in this DG0 case. You’ll get suboptimal convergence with the P2P0 element anyway, but it’ll be interesting to see how @dokken’s suggestion performs.

1 Like
4

Thank you for your answer @dokken, I had forgotten that the degree of freedom for a DG function was at the cell center. If the unit square is subdivided into m subdivisions in each direction, a cell will be a rectangular isosceles triangle of size 1/m. So the distance between a vertex of a cell and its center should be less than 0.5/m.

Following this idea, I tried to set p(0.0, 0.0) = 0 (which is equivalent with setting p to 0 on the whole boundary in my case) in the code below, but unfortunately I get the same problem :neutral_face:

from dolfin import *

# Load mesh 
mesh_subdiv = 32
mesh = UnitSquareMesh(mesh_subdiv, mesh_subdiv)

# Define function spaces
Elem_u = VectorElement("CG", mesh.ufl_cell(), 2)
Elem_p = FiniteElement("DG", mesh.ufl_cell(), 0)
W = FunctionSpace(mesh, Elem_u * Elem_p)

# Define Dirichlet boundary
def velocity_boundary(x):
    return x[0] < DOLFIN_EPS or x[0] > 1.0 - DOLFIN_EPS or x[1] < DOLFIN_EPS or x[1] > 1.0 - DOLFIN_EPS

tol = 0.5/mesh_subdiv
def pressure_boundary(x, tol):
    return abs(x[0]) < tol and abs(x[1]) < tol

# Boundary conditions for velocity and pressure
bc0 = DirichletBC(W.sub(0), Constant((0.0, 0.0)), velocity_boundary)        
bc1 = DirichletBC(W.sub(1), Constant(0.0), pressure_boundary, "pointwise")
bcs = [bc0, bc1]

# Define variational problem
(u, p) = TrialFunctions(W)
(v, q) = TestFunctions(W)
f = Expression(("-4*pi*pi*sin(2*pi*x[1])","-16*pi*pi*sin(2*pi*x[0])*cos(2*pi*x[1]) + 4*pi*pi*sin(2*pi*x[0])"), degree=2)
a = (inner(grad(u), grad(v)) - div(v)*p + q*div(u))*dx
L = inner(f, v)*dx

# Compute solution
w = Function(W)
solve(a == L, w, bcs)

@nate what do you mean by “weak imposition of the BCs” ? I would prefer to avoid adding the constraint \int_{\Omega} p \ \mathrm{d}x = 0 if possible.

5

Hi,

Here is a little update: I solved my problem by removing the boundary condition on the pressure and adding a penalty-term \varepsilon \int_{\Omega} pq in the variational formulation. For \varepsilon small enough, this gives the correct pressure :slight_smile:

Yet, very small \varepsilon yields to a badly conditioned assembled matrix, so I would be curious to know whether there exists a solution using only DirichletBC and not changing the formulation.