# Steady Stokes equation (3D) - dolfinx

**URL:** <https://fenicsproject.discourse.group/t/steady-stokes-equation-3d-dolfinx/9709>\
**Category:** General\
**Created:** [November 8, 2022, 4:46pm UTC](https://fenicsproject.discourse.group/t/steady-stokes-equation-3d-dolfinx/9709 "2022-11-08T16:46:32Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![Stefano\_Riva](https://avatars.discourse-cdn.com/v4/letter/s/cc9497/32.png) [@Stefano\_Riva](https://fenicsproject.discourse.group/u/Stefano_Riva)\
**Post date:** [November 8, 2022, 4:46pm UTC](https://fenicsproject.discourse.group/t/steady-stokes-equation-3d-dolfinx/9709/1 "2022-11-08T16:46:32Z")

</div>

Hi everyone,  
I’m using FEniCSx to implement a solver for 3D Stokes equation

\nu\Delta \mathbf{u}-\nabla p = 0

I have implemented some lines of code to deal with the 2D version of the problem and the results are reasonable. However, when moving to 3D nothing seems to work, there is no error in output but the solution is a nonsense … I have modified the mesh generation and the assignment of the BC.

Do you have some hints or tutorials in which a similar problem has been studied (also steady Navier-Stokes would be helpful)?

Thanks

PS: I’m attaching the link to the [2D code](https://github.com/Steriva/FEniCSx/blob/main/2D/BFSsym_2D.ipynb).

---

<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:** [November 8, 2022, 4:57pm UTC](https://fenicsproject.discourse.group/t/steady-stokes-equation-3d-dolfinx/9709/2 "2022-11-08T16:57:40Z")

</div>

You are stating that you want to solve the Stokes equations, which for instance has been handled in [Stokes equations with Taylor-Hood elements — DOLFINx 0.5.1 documentation](https://docs.fenicsproject.org/dolfinx/v0.5.1/python/demos/demo_stokes.html)

---

<div class="post-metadata">

**Author:** ![Stefano\_Riva](https://avatars.discourse-cdn.com/v4/letter/s/cc9497/32.png) [@Stefano\_Riva](https://fenicsproject.discourse.group/u/Stefano_Riva)\
**Post date:** [November 8, 2022, 5:02pm UTC](https://fenicsproject.discourse.group/t/steady-stokes-equation-3d-dolfinx/9709/3 "2022-11-08T17:02:45Z")

</div>

Yes, I have looked at that demo and I have implemented the attached code for 2D problem. However, I’m not able to extent to 3D

---

<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:** [November 8, 2022, 9:16pm UTC](https://fenicsproject.discourse.group/t/steady-stokes-equation-3d-dolfinx/9709/4 "2022-11-08T21:16:39Z")

</div>

It is quite straightforward to port demos from DOLFIN to DOLFINx.  
For instance, take:  
[https://fenicsproject.org/olddocs/dolfin/latest/python/demos/stokes-iterative/demo\_stokes-iterative.py.html](https://fenicsproject.org/olddocs/dolfin/latest/python/demos/stokes-iterative/demo_stokes-iterative.py.html)  
This can be converted as follows:

```python
# Demo solving Stokes problem in 3D in DOLFINx using MINRES and AMG
# Copyright: 2022
# Author: Jørgen S. Dokken
# License: MIT

import numpy as np
import dolfinx
import ufl
from petsc4py import PETSc
from mpi4py import MPI

N = 16
mesh = dolfinx.mesh.create_unit_cube(MPI.COMM_WORLD, N, N, N,
                                     cell_type=dolfinx.mesh.CellType.hexahedron)
P2 = ufl.VectorElement("Lagrange", mesh.ufl_cell(), 2)
P1 = ufl.FiniteElement("Lagrange", mesh.ufl_cell(), 1)
TH = ufl.MixedElement([P2, P1])
W = dolfinx.fem.FunctionSpace(mesh, TH)

def right(x, tol=1e-14): return x[0] > 1 - tol

def top_bottom(x, tol=1e-14):
    return np.logical_or(
        np.isclose(x[1], 0, atol=tol), np.isclose(x[1], 1, atol=tol))

V, V_to_W = W.sub(0).collapse()

# No-slip boundary condition for velocity
nonslip = dolfinx.fem.Function(V)
nonslip.x.set(0.0)
nonslip_facets = dolfinx.mesh.locate_entities_boundary(
    mesh, mesh.topology.dim-1, top_bottom)
nonslip_dofs = dolfinx.fem.locate_dofs_topological(
    (W.sub(0), V), mesh.topology.dim-1, nonslip_facets)
bc0 = dolfinx.fem.dirichletbc(nonslip, nonslip_dofs, W.sub(0))

# Inflow boundary condition for velocity

def inflow(x):
    zero_comp = np.zeros(x.shape[1], dtype=np.float64)
    return (-np.sin(np.pi*x[1]), zero_comp, zero_comp)

inflow_facets = dolfinx.mesh.locate_entities_boundary(mesh, mesh.topology.dim-1,
                                                      right)
inflow_dofs = dolfinx.fem.locate_dofs_topological((W.sub(0), V), mesh.topology.dim-1,
                                                  inflow_facets)
u_bc = dolfinx.fem.Function(V)
u_bc.interpolate(inflow)
bc1 = dolfinx.fem.dirichletbc(u_bc, inflow_dofs, W.sub(0))

bcs = [bc0, bc1]

(u, p) = ufl.TrialFunctions(W)
(v, q) = ufl.TestFunctions(W)
f = dolfinx.fem.Constant(mesh, (0., 0., 0.))

dx = ufl.Measure("dx", domain=mesh)
a = ufl.inner(ufl.grad(u), ufl.grad(v)) * dx + \
    ufl.div(v)*p*dx + q*ufl.div(u)*dx
L = ufl.inner(f, v) * dx

a_f = dolfinx.fem.form(a)
L_f = dolfinx.fem.form(L)
A = dolfinx.fem.petsc.assemble_matrix(a_f, bcs=bcs)
A.assemble()
b = dolfinx.fem.petsc.assemble_vector(L_f)
dolfinx.fem.petsc.apply_lifting(b, [a_f], [bcs])
b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
dolfinx.fem.petsc.set_bc(b, bcs)
b.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)

# Preconditioner matrix form
Pf = dolfinx.fem.form(ufl.inner(ufl.grad(u), ufl.grad(v)) * dx + p * q * dx)
P = dolfinx.fem.petsc.assemble_matrix(Pf, bcs=bcs)
P.assemble()

ksp = PETSc.KSP().create(mesh.comm)
ksp.setOperators(A, P)
ksp.setType(PETSc.KSP.Type.MINRES)
pc = ksp.getPC()
pc.setType("hypre")
pc.setHYPREType("boomeramg")
pc.setUp()

# To view Preconditioner info
# pc.view()

wh = dolfinx.fem.Function(W)
ksp.solve(b, wh.vector)

# To view solver info
# ksp.view()

uh = wh.sub(0).collapse()
p = wh.sub(1).collapse()

print(f"Converged Reason {ksp.getConvergedReason()}" +
      f"\nNumber of iterations {ksp.getIterationNumber()}")

with dolfinx.io.VTXWriter(mesh.comm, "u.bp", [uh]) as vtx:
    vtx.write(0.)

with dolfinx.io.VTXWriter(mesh.comm, "p.bp", [uh]) as vtx:
    vtx.write(0.)

```

---

<div class="post-metadata">

**Author:** ![Stefano\_Riva](https://avatars.discourse-cdn.com/v4/letter/s/cc9497/32.png) [@Stefano\_Riva](https://fenicsproject.discourse.group/u/Stefano_Riva)\
**Post date:** [November 9, 2022, 9:51am UTC](https://fenicsproject.discourse.group/t/steady-stokes-equation-3d-dolfinx/9709/5 "2022-11-09T09:51:10Z")

</div>

Thank you very much for the support, Mr. Dokken. I have translated the code you gave me into my case and it works well.

I have one additional question: is it possible to use a P1b-P1 couple of functional space, in order to reduce the dimension of the system and to avoid the implementation of stabilisation techniques?  
I have found a discussion about this for scalar function, does it works the same for vector quantities?

---

<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:** [November 9, 2022, 10:08am UTC](https://fenicsproject.discourse.group/t/steady-stokes-equation-3d-dolfinx/9709/6 "2022-11-09T10:08:26Z")

</div>

See for instance: [https://kent-and.github.io/mek4250/lecture\_08\_stokes\_alt.pdf](https://kent-and.github.io/mek4250/lecture_08_stokes_alt.pdf) (slide 9) and [The Stokes equations — FEniCS 22 tutorial](https://jorgensd.github.io/fenics22-tutorial/comparing_elements.html#piecewise-linear-pressure-space)