# Convergence in Stokes problem

**URL:** <https://fenicsproject.discourse.group/t/convergence-in-stokes-problem/11924>\
**Category:** variational formulation\
**Created:** [August 8, 2023, 2:48pm UTC](https://fenicsproject.discourse.group/t/convergence-in-stokes-problem/11924 "2023-08-08T14:48:54Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![michele\_barucca](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/michele_barucca/32/5350_2.png) [@michele\_barucca](https://fenicsproject.discourse.group/u/michele_barucca)\
**Post date:** [August 8, 2023, 2:48pm UTC](https://fenicsproject.discourse.group/t/convergence-in-stokes-problem/11924/1 "2023-08-08T14:48:54Z")

</div>

Good evening,

in the Stokes problem using P1-P1 finite elements with a stabilization term (I implemented the easiest, i.e. the Brezzi-Pitkaranta stabilization), one should expect a convergence like O(h) for the H1 error in the velocity and the pressure. In my case this does not happens and following previous advice I tried first to solve this problem using a direct solver, but still the results are not right.  
The minimal code is the following:

```auto
import numpy as np

import ufl
from dolfinx import cpp as _cpp
from dolfinx import fem
from dolfinx.fem import (Constant, Function, FunctionSpace, dirichletbc,
                         extract_function_spaces, form,
                         locate_dofs_geometrical, locate_dofs_topological)
from dolfinx.io import XDMFFile
from dolfinx.mesh import (CellType, GhostMode, create_rectangle, create_unit_cube,
                          locate_entities_boundary, exterior_facet_indices)
from ufl import div, dx, grad, inner

from mpi4py import MPI
from petsc4py import PETSc

# Create mesh
N = 32
msh = create_unit_cube(MPI.COMM_WORLD,
                       N, N, N,
                       CellType.hexahedron, GhostMode.none)

'''# Function to mark x = 0, x = 1 and y = 0
def noslip_boundary(x):
    return np.logical_or(np.logical_or(np.isclose(x[0], 0.0),
                                       np.isclose(x[0], 1.0)),
                         np.isclose(x[1], 0.0))

# Function to mark the lid (y = 1)
def lid(x):
    return np.isclose(x[1], 1.0)

# Lid velocity
def lid_velocity_expression(x):
    return np.stack((np.ones(x.shape[1]), np.zeros(x.shape[1])))'''

def u_ex(x): 
        return np.vstack((4*(x[1]-1)*x[1]*(x[0]-1)*(x[0]-1)*x[0]*x[0]*(x[1]-0.5), -4*(x[0]-0.5)*x[1]*x[1]*x[0]*(x[1]-1)*(x[1]-1)*(x[0]-1), np.zeros_like(x[0])))
def p_ex(x):
    return x[0]+x[1]+x[2]-1.5

def f_ex(x): 
    return np.vstack((1-24*(-1 + x[0])*(x[0]-1)*x[0]*x[0]*(-0.5 + x[1])- 48*(0.166667 - x[0] + x[0]*x[0])*(-1 + x[1])*(-0.5 + x[1])*x[1], 1+24*(-0.5 + x[0])*(-1 + x[1])*(x[1]-1)*x[1]*x[1] + 8*(-1 + x[0])*(-0.5 + x[0])*x[0]*(1 - 6*x[1] + 6*x[1]*x[1]), np.ones_like(x[0])))

P2 = ufl.VectorElement("Lagrange", msh.ufl_cell(), 1)
P1 = ufl.FiniteElement("Lagrange", msh.ufl_cell(), 1)

# Create the function space
TH = P2 * P1
W = FunctionSpace(msh, TH)
W0, _ = W.sub(0).collapse()
W1, _ = W.sub(1).collapse()

'''# No slip boundary condition
noslip = Function(V)
facets = locate_entities_boundary(msh, 1, noslip_boundary)
dofs = locate_dofs_topological((W.sub(0), V), 1, facets)
bc0 = dirichletbc(noslip, dofs, W.sub(0))

# Driving velocity condition u = (1, 0) on top boundary (y = 1)
lid_velocity = Function(W0)
lid_velocity.interpolate(lid_velocity_expression)
facets = locate_entities_boundary(msh, 1, lid)
dofs = locate_dofs_topological((W.sub(0), V), 1, facets)
bc1 = dirichletbc(lid_velocity, dofs, W.sub(0))'''

uD = Function(W0)
uD.interpolate(u_ex)
uD.vector.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
# Create facet to cell connectivity required to determine boundary facets
tdim = msh.topology.dim
fdim = tdim - 1
msh.topology.create_connectivity(fdim, tdim)
boundary_facets = exterior_facet_indices(msh.topology)
boundary_dofs = locate_dofs_topological((W.sub(0),W0), fdim, boundary_facets)
bc1 = dirichletbc(uD, boundary_dofs, W.sub(0))

# Since for this problem the pressure is only determined up to a
# constant, we pin the pressure at the point (0, 0)
zero = Function(W1)
zero.x.set(-1.5)
zero.vector.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
dofs = locate_dofs_geometrical((W.sub(1), W1), lambda x: np.isclose(x.T, [0, 0, 0]).all(axis=1))
bc2 = dirichletbc(zero, dofs, W.sub(1))

# Collect Dirichlet boundary conditions
bcs = [bc1, bc2]

# Define variational problem
(u, p) = ufl.TrialFunctions(W)
(v, q) = ufl.TestFunctions(W)
f = Function(W0)
f.interpolate(f_ex)
h = PETSc.ScalarType(ufl.sqrt(2)*1./N)
beta = PETSc.ScalarType(0.1)
a = form((inner(grad(u), grad(v)) + inner(p, div(v)) + inner(div(u), q)) * dx - h*h*beta*ufl.inner(ufl.nabla_grad(p),ufl.nabla_grad(q))*dx)
L = form(inner(f, v) * dx)

# Assemble LHS matrix and RHS vector
A = fem.petsc.assemble_matrix(a, bcs=bcs)
A.assemble()
b = fem.petsc.assemble_vector(L)

fem.petsc.apply_lifting(b, [a], bcs=[bcs])
b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)

# Set Dirichlet boundary condition values in the RHS
fem.petsc.set_bc(b, bcs)

# Create and configure solver
ksp = PETSc.KSP().create(msh.comm)
ksp.setOperators(A)
ksp.setType("preonly")
ksp.getPC().setType("lu")
ksp.getPC().setFactorSolverType("mumps")
opts = PETSc.Options()
opts["mat_mumps_icntl_14"] = 80 # Increase MUMPS working memory
opts["mat_mumps_icntl_24"] = 0 # Option to support solving a singular matrix (pressure nullspace)
opts["mat_mumps_icntl_25"] = 0 # Option to support solving a singular matrix (pressure nullspace)
opts["ksp_error_if_not_converged"] = 1
opts["ksp_monitor"] = None
ksp.setFromOptions()

# Compute the solution
U = Function(W)
ksp.solve(b, U.vector)

# Split the mixed solution and collapse
u = U.sub(0).collapse()
p = U.sub(1).collapse()

L2_error = fem.form(ufl.inner(u - uD, u - uD) * ufl.dx)
error_local = fem.assemble_scalar(L2_error)
error_L2 = np.sqrt(msh.comm.allreduce(error_local, op=MPI.SUM))

if msh.comm.rank == 0:
    print(f"Error_L2 : {error_L2:.2e}")

eh = u - uD
error_H10 = fem.form(ufl.inner(ufl.grad(eh), ufl.grad(eh)) * ufl.dx)
E_H10 = np.sqrt(msh.comm.allreduce(fem.assemble_scalar(error_H10), op=MPI.SUM))
if msh.comm.rank == 0:
    print(f"H01-error: {E_H10:.2e}")

pD = Function(W1)
pD.interpolate(p_ex)
L2_errorp = fem.form(ufl.inner(p - pD, p - pD) * ufl.dx)
error_localp = fem.assemble_scalar(L2_errorp)
error_L2p = np.sqrt(msh.comm.allreduce(error_localp, op=MPI.SUM))

if msh.comm.rank == 0:
    print(f"Error_L2p : {error_L2p:.2e}")

# Write the solution to file
with XDMFFile(MPI.COMM_WORLD, "out_stokes/u.xdmf", "w") as ufile_xdmf:
    ufile_xdmf.write_mesh(msh)
    ufile_xdmf.write_function(u)

with XDMFFile(MPI.COMM_WORLD, "out_stokes/p.xdmf", "w") as pfile_xdmf:
    pfile_xdmf.write_mesh(msh)
    pfile_xdmf.write_function(p)

with XDMFFile(MPI.COMM_WORLD, "out_stokes/uex.xdmf", "w") as pfile_xdmf:
    pfile_xdmf.write_mesh(msh)
    pfile_xdmf.write_function(uD)

```

The error for a mesh with respectively 8, 16, 32 elements are:

```auto
Error_L2 : 5.08e-03
H01-error: 3.39e-02
Error_L2p : 2.67e+00

Error_L2 : 1.41e-03
H01-error: 1.14e-02
Error_L2p : 2.90e+00

Error_L2 : 3.68e-04
H01-error: 3.71e-03
Error_L2p : 3.03e+00

```

Which are note corrects (in particular the one of the pressure??). Can someone give me an hand on this problem that is driving me crazy?

---

<div class="post-metadata">

**Author:** ![nate](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/nate/32/17_2.png) [@nate](https://fenicsproject.discourse.group/u/nate)\
**Post date:** [August 8, 2023, 3:02pm UTC](https://fenicsproject.discourse.group/t/convergence-in-stokes-problem/11924/2 "2023-08-08T15:02:03Z")

</div>

Looks like you have the wrong sign in front of the pressure term.

```python
a = form((inner(grad(u), grad(v)) + inner(p, div(v)) + inner(div(u), q)) * dx - h*h*beta*ufl.inner(ufl.nabla_grad(p),ufl.nabla_grad(q))*dx)

```

You could either flip the signs on the continuity and pressure terms to maintain symmetry, or just post process your pressure solution to account for it.

---

<div class="post-metadata">

**Author:** ![michele\_barucca](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/michele_barucca/32/5350_2.png) [@michele\_barucca](https://fenicsproject.discourse.group/u/michele_barucca)\
**Post date:** [August 8, 2023, 3:11pm UTC](https://fenicsproject.discourse.group/t/convergence-in-stokes-problem/11924/3 "2023-08-08T15:11:33Z")

</div>

Thanks for the answer. This solved the bed convergence of the pressure, but still I don’t get the rate of convergence I expect.  
Did you notice something else that is wrong?