# Singular Stokes Problem

**URL:** <https://fenicsproject.discourse.group/t/singular-stokes-problem/11052>\
**Category:** General\
**Created:** [April 28, 2023, 2:22pm UTC](https://fenicsproject.discourse.group/t/singular-stokes-problem/11052 "2023-04-28T14:22:13Z")\
**Posts on this page:** 7\
**Page:** 1

<div class="post-metadata">

**Author:** ![ahandeep\_manna](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/ahandeep_manna/32/2897_2.png) [@ahandeep\_manna](https://fenicsproject.discourse.group/u/ahandeep_manna)\
**Post date:** [April 28, 2023, 2:22pm UTC](https://fenicsproject.discourse.group/t/singular-stokes-problem/11052/1 "2023-04-28T14:22:13Z")

</div>

Dear All,

I’m trying to solve the stokes equation with Pure Neumann boundary conditions or Traction boundary conditions at the boundaries.

 ![Screen Shot 2023-04-28 at 16.18.51](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/2/29cca361fe859c9af82f0a1d16b2818666f22072.png)  
This problem is singular and similiar to ([https://fenicsproject.org/olddocs/dolfin/1.3.0/python/demo/documented/singular-poisson/python/documentation.html#](https://fenicsproject.org/olddocs/dolfin/1.3.0/python/demo/documented/singular-poisson/python/documentation.html#)) because the velocity field and the pressure field are not determined up to a constant .

> mesh=UnitSquareMesh(20,20)  
> V = VectorFunctionSpace(mesh,‘P’, 1)  
> p\_element=FiniteElement(‘P’, mesh.ufl\_cell(),1)  
> v\_element=VectorElement(‘P’,mesh.ufl\_cell(),2)  
> element = MixedElement([p\_element, v\_element])  
> FS=FunctionSpace(mesh,element)
> 
> q,v=TestFunctions(FS)  
> P,u=TrialFunctions(FS)
> 
> force=df.Expression((‘x[0]’,‘0’),degree=1)  
> a = inner(grad(u), grad(v))_df.dx + div(v)pdf.dx + q_div(u)\*df.dx  
> L=inner(force,v)\*df.dx

1. Can someone help me with how to remove the null space of the Matrix (a), such that i can eliminate constant velocity and pressure solutions ? similiar to ([https://fenicsproject.org/qa/4121/help-parallelising-space-method-eliminating-rigid-motion/](https://fenicsproject.org/qa/4121/help-parallelising-space-method-eliminating-rigid-motion/)) .
2. For the elasticity problem however there is only one function space. How can one do this procedure for both the pressure and the velocity?

---

<div class="post-metadata">

**Author:** ![Kei\_Yamamoto](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/kei_yamamoto/32/2904_2.png) [@Kei\_Yamamoto](https://fenicsproject.discourse.group/u/Kei_Yamamoto)\
**Post date:** [May 1, 2023, 7:50am UTC](https://fenicsproject.discourse.group/t/singular-stokes-problem/11052/2 "2023-05-01T07:50:10Z")

</div>

Have you looked at this post ?  
[https://fenicsproject.org/qa/13927/pressure-nullspace-for-stokes-problem/](https://fenicsproject.org/qa/13927/pressure-nullspace-for-stokes-problem/)

---

<div class="post-metadata">

**Author:** ![JulietteD](https://avatars.discourse-cdn.com/v4/letter/j/a9adbd/32.png) [@JulietteD](https://fenicsproject.discourse.group/u/JulietteD)\
**Post date:** [March 25, 2025, 2:50pm UTC](https://fenicsproject.discourse.group/t/singular-stokes-problem/11052/3 "2025-03-25T14:50:41Z")

</div>

Hello,

I am also trying to solve the Stokes equation with pure Neumann boundary condition, but I don’t manage to define the correct nullspace.  
The problem I am trying to solve is

\begin{align} \Delta {\bf u} - \nabla p = 0 \text{ in } \Omega \\ \text{div} \, {\bf u} = 0 \text{ in } \Omega \\ (\nabla {\bf u}) \cdot {\bf n} - p \, {\bf n} = {\bf g} \text{ on } \partial \Omega \end{align}

and g satisfies the compatibility condition \int\_{\partial \Omega} {\bf g} = 0.

Here is the minimal working example with a tentative nullspace:

```auto
from mpi4py import MPI
from petsc4py import PETSc
dtype = PETSc.ScalarType  
import ufl
from basix.ufl import element, mixed_element
import dolfinx
from dolfinx.fem.petsc import assemble_matrix_block, assemble_vector_block

domain = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 40, 40)

V_el = element("Lagrange", domain.basix_cell(), 2, shape=(2,)) 
Q_el = element("Lagrange", domain.basix_cell(), 1) 
W_el = mixed_element([V_el,Q_el])
W = dolfinx.fem.functionspace(domain, W_el)
V,V_to_W = W.sub(0).collapse()
Q,Q_to_W = W.sub(1).collapse()

## Definition of the linear forms
(u, p) = ufl.TrialFunction(V), ufl.TrialFunction(Q)
(v, q) = ufl.TestFunction(V), ufl.TestFunction(Q)
a = dolfinx.fem.form([[ufl.inner(ufl.grad(u), ufl.grad(v)) * ufl.dx , ufl.inner(p, ufl.div(v)) * ufl.dx], [ufl.inner(ufl.div(u), q) * ufl.dx, None]])

## Neumann BC
x = ufl.SpatialCoordinate(domain)
g = ufl.as_vector([2*x[0]-1, 0.2*x[1]-0.1])
L = dolfinx.fem.form([ufl.inner(g, v) * ufl.ds, ufl.inner(dolfinx.fem.Constant(domain, dtype(0)), q) * ufl.dx])

# Assemble the block operator and RHS vector
A = assemble_matrix_block(a)
A.assemble()
b = assemble_vector_block(L, a)

## Set the nullspace 
bs = W.dofmap.index_map_bs
length0 = W.dofmap.index_map.size_local

basis = [dolfinx.la.vector(W.dofmap.index_map, bs=bs, dtype=dtype) for i in range(3)]
basis_arrays = [_b.array for _b in basis]
dofs_v = [V.sub(i).dofmap.list.flatten() for i in range(2)]
dofs_p = [Q.dofmap.list.flatten()]
basis_arrays[0][dofs_v[0]] = 1.0
basis_arrays[1][dofs_v[1]] = 1.0
basis_arrays[2][dofs_p] = 1.0
dolfinx.la.orthonormalize(basis)
basis_petsc = [
    PETSc.Vec().createWithArray(x[: bs * length0], bsize=bs, comm=W.mesh.comm) 
    for x in basis_arrays
]
nullspace = PETSc.NullSpace().create(vectors=basis_petsc)  

assert nullspace.test(A)
A.setNearNullSpace(nullspace)

# Create a block vector (x) to store the full solution, and solve
# Create a solver
ksp = PETSc.KSP().create(domain.comm)
ksp.setOperators(A)
ksp.setType("preonly")

# Set the solver type to MUMPS (LU solver)
pc = ksp.getPC()
pc.setType("lu")
pc.setFactorSolverType("superlu_dist")
pc.setFactorSetUpSolverType()

# Create a block vector (x) to store the full solution, and solve
x = A.createVecLeft()
ksp.solve(b, x)

# Create Functions and scatter x solution
u, p = dolfinx.fem.Function(V), dolfinx.fem.Function(Q)
offset = V.dofmap.index_map.size_local * V.dofmap.index_map_bs
u.x.array[:offset] = x.array_r[:offset]
p.x.array[: (len(x.array_r) - offset)] = x.array_r[offset:]

```

The `assert nullspace.test(A)` raises an error. Is there something wrong with the way I am creating the nullspace? I am quite new to dolfinx.

Any help is appreciated!

---

<div class="post-metadata">

**Author:** ![Stein](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/stein/32/8751_2.png) [@Stein](https://fenicsproject.discourse.group/u/Stein)\
**Post date:** [March 25, 2025, 3:44pm UTC](https://fenicsproject.discourse.group/t/singular-stokes-problem/11052/4 "2025-03-25T15:44:04Z")

</div>

In your problem formulation, there is no remaining null space. The constant of the pressure is fixed due to the natural boundary condition: p+const no-longer satisfies the PDE.

---

<div class="post-metadata">

**Author:** ![JulietteD](https://avatars.discourse-cdn.com/v4/letter/j/a9adbd/32.png) [@JulietteD](https://fenicsproject.discourse.group/u/JulietteD)\
**Post date:** [March 25, 2025, 4:14pm UTC](https://fenicsproject.discourse.group/t/singular-stokes-problem/11052/5 "2025-03-25T16:14:30Z")

</div>

Thank you for your quick answer! When you say that there is no nullspace, do you mean only for the pressure or also for the velocity ? The problem seems to be satisfied by u and by u + cst.

I tried changing the nullspace to only take into account constant velocities, i.e. keep the first two `basis` vectors. The nullspace `assert` does not raise any errors anymore, but the solution does not converge when I refine the mesh. Is this a well-known issue and if so, do you have any pointers?

---

<div class="post-metadata">

**Author:** ![Stein](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/stein/32/8751_2.png) [@Stein](https://fenicsproject.discourse.group/u/Stein)\
**Post date:** [March 25, 2025, 4:51pm UTC](https://fenicsproject.discourse.group/t/singular-stokes-problem/11052/6 "2025-03-25T16:51:49Z")

</div>

Ah, you’re right. For most flow problems (with advection operators and inflow conditions) it is the pressure nullspace that is the issue, hence my automatic fixation. This relates more closely to the null-spaces observed in solids. Is this demo of any use to you? [Elasticity using algebraic multigrid — DOLFINx 0.9.0 documentation](https://docs.fenicsproject.org/dolfinx/v0.9.0/python/demos/demo_elasticity.html)

---

<div class="post-metadata">

**Author:** ![JulietteD](https://avatars.discourse-cdn.com/v4/letter/j/a9adbd/32.png) [@JulietteD](https://fenicsproject.discourse.group/u/JulietteD)\
**Post date:** [March 25, 2025, 6:33pm UTC](https://fenicsproject.discourse.group/t/singular-stokes-problem/11052/7 "2025-03-25T18:33:48Z")

</div>

Thank you for the pointer. I tried adapting the example to the 2D case:

```auto
## Set the nullspace 
bs = W.dofmap.index_map_bs
length0 = W.dofmap.index_map.size_local

basis = [dolfinx.la.vector(W.dofmap.index_map, bs=bs, dtype=dtype) for i in range(3)]
basis_arrays = [_b.array for _b in basis]
dofs_v = [V.sub(i).dofmap.list.flatten() for i in range(2)]
# translation modes
basis_arrays[0][dofs_v[0]] = 1.0
basis_arrays[1][dofs_v[1]] = 1.0

# rotation modes
xdof = V.tabulate_dof_coordinates()
dofs_block = V.dofmap.list.flatten()
x0, x1= xdof[dofs_block, 0], xdof[dofs_block, 1]
basis_arrays[2][dofs_v[0]] = - x1
basis_arrays[2][dofs_v[1]] = x0 

dolfinx.la.orthonormalize(basis)
basis_petsc = [
    PETSc.Vec().createWithArray(_x[: bs * length0], bsize=bs, comm=W.mesh.comm) 
    for _x in basis_arrays
]
nullspace = PETSc.NullSpace().create(vectors=basis_petsc)  

```

But the assert raises an error again.  
From a mathematical point of view though, including the rigid body motions in the nullspace seems to be the good approach. Maybe I am not defining them correctly?