Stokes equation with pressure dirichlet condition

Hello, I want to solve Stoke equation with pressure dirichlet condition.

The geometry of my case is as followed.

The left side is inlet (denoted by blue) and the right side is outlet (denoted by green), and the rest boundaries is noslip wall (denoted by red).

The inlet dirichlet condition is p|_{left}=2, and the outlet dirichlet condition is p|_{right}=0.

And I refer to Stokes equations using Taylor-Hood elements , the code is followed.

import dolfinx.fem.petsc
import gmsh
import numpy as np
import ufl
import mpi4py
from basix.ufl import element, mixed_element
from dolfinx import fem
from dolfinx.fem import Function, dirichletbc, form, functionspace, locate_dofs_topological
from ufl import div, dx, grad, inner, ds, dot, nabla_grad, sym, Identity, FacetNormal
import viskex
import dolfinx
from petsc4py import PETSc

# geometry & mesh parameter (unit is m)
L, meshsz = 0.001, 2E-5
WIDTH, HEIGHT = 5 * L, L

# fluid parameter
rho, mu = 1E3, 1E-3  # kg/m^3, Pa * s
init_pressure = 2  # Pa

gmsh.initialize()
gmsh.option.setNumber("Mesh.MshFileVersion", 2.2)
gmsh.option.setNumber("Mesh.Algorithm", 5)
gmsh.option.setNumber('General.Terminal', 0)
gmsh.option.setNumber('Mesh.MeshSizeMax', meshsz)

geo = gmsh.model.occ.addRectangle(-WIDTH/2, 0, 0, WIDTH*2, HEIGHT)
gmsh.model.occ.synchronize()
gmsh.model.mesh.generate(2)

bd = gmsh.model.occ.getEntities(1)
wall, inlet, outlet = [], [], []
for reg in bd:
    tag = reg[1]
    x, y, _ = gmsh.model.occ.getCenterOfMass(1, tag)
    if x < -WIDTH / 2 + meshsz / 2:
        inlet.append(tag)
    elif x > 3 * WIDTH / 2 - meshsz / 2:
        outlet.append(tag)
    else:
        wall.append(tag)
gmsh.model.addPhysicalGroup(1, inlet, 1, name="inlet")
gmsh.model.addPhysicalGroup(1, outlet, 2, name="outlet")
gmsh.model.addPhysicalGroup(1, wall, 3, name="wall")
gmsh.model.addPhysicalGroup(2, [geo], 0, name="solid")

mesh, subdomains, boundaries = dolfinx.io.gmshio.model_to_mesh(gmsh.model, comm=mpi4py.MPI.COMM_WORLD, rank=0, gdim=2)
mesh.topology.create_connectivity(mesh.topology.dim - 1, mesh.topology.dim)
inlet = boundaries.indices[boundaries.values == 1]
outlet = boundaries.indices[boundaries.values == 2]
wall = boundaries.indices[boundaries.values == 3]
viskex.dolfinx.plot_mesh_tags(mesh, boundaries, "boundaries")

# refer https://docs.fenicsproject.org/dolfinx/main/python/demos/demo_stokes.html
def mixed_direct(msh, inlet, outlet, wall):

    def inlet_pressure_expression(x):
        return np.ones(x.shape[1]) * init_pressure

    P2 = element("Lagrange", msh.basix_cell(), 2, shape=(msh.geometry.dim,))
    P1 = element("Lagrange", msh.basix_cell(), 1)

    TH = mixed_element([P2, P1])
    W = functionspace(msh, TH)
    W0, W1 = W.sub(0), W.sub(1)
    V, _ = W0.collapse()
    Q, _ = W1.collapse()

    # No slip boundary condition
    noslip = Function(V)
    dofs = locate_dofs_topological((W0, V), 1, wall)
    bc0 = dirichletbc(noslip, dofs, W0)

    # inlet boundary condition: pressure = 2Pa
    inlet_pressure = Function(Q)
    inlet_pressure.interpolate(inlet_pressure_expression)
    dofs = locate_dofs_topological((W1, Q), 1, inlet)
    bc1 = dirichletbc(inlet_pressure, dofs, W1)

    # outlet boundary condition: pressure = 0Pa
    outlet_pressure = Function(Q)
    dofs = locate_dofs_topological((W1, Q), 1, outlet)
    bc2 = dirichletbc(outlet_pressure, dofs, W1)

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

    # Define variational problem
    (u, p) = ufl.TrialFunctions(W)
    (v, q) = ufl.TestFunctions(W)
    f = Function(V)
    a = form((mu * inner(grad(u), grad(v)) + inner(p, div(v)) + inner(div(u), q)) * dx)
    L = form(inner(f, v) * dx)

    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")

    # Configure MUMPS to handle pressure nullspace
    pc = ksp.getPC()
    pc.setType("lu")
    pc.setFactorSolverType("mumps")
    pc.setFactorSetUpSolverType()
    pc.getFactorMatrix().setMumpsIcntl(icntl=24, ival=1)
    pc.getFactorMatrix().setMumpsIcntl(icntl=25, ival=0)

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

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

    viskex.dolfinx.plot_vector_field(u, "u", glyph_factor=1e-2)
    viskex.dolfinx.plot_scalar_field(p, "p")


mixed_direct(mesh, inlet=inlet, outlet=outlet, wall=wall)

But the result seems to be not right, the velocity field is shown as followed

The velocity seem to have value only at the inlet boundary, the rest is almost 0.

I think there should be error in my code, but I’m not familiar to dolfinx.

The version of dolfinx I used is dolfinx0.8.0.

There might be an error in the code, but typically in FEM we never impose a Dirichlet BC on the pressure. If you write down the weak formulation on paper you’ll realize that the pressure isn’t derived any more, which gives you no room to impose a Dirichlet BC on that variable (since there are no more differential operators on it).

I have not run the code itself, but if there was an error there I would highlight that you seem to be misusing mumps options. You have a part with Configure MUMPS to handle pressure nullspace which at the very least is definitely not needed, since the pressure does not have a null space with this BCs, and worst case scenario might be detrimental to the solver.

Sorry, I don’t have time for a complete answer, but it is possible to impose a condition on the pressure for both inlet and outlet and get a well defined problem.

For the Dirichlet condition you requested, take a look at Theory and Practice of Finite Elements | SpringerLink page 182 - there valid boundary conditions based on the curl-curl formulation of the vector Laplacian are derived. However, you will get the conditions

u \times n = 0 \text{ on } \Gamma \quad \text{ and } \quad (p - \mu(u\cdot n)) = -g \text{ on } \Gamma

for some \mu_3 \geq 0. So I guess this is what you really want as boundary condition. The second BC which involves the pressure can then be implemented weakly.

Another possibility is to use the formulation considered in https://doi.org/10.1002/(SICI)1097-0363(19960315)22:5<325::AID-FLD307>3.0.CO;2-Y
There, it is shown how you can “modify” the standard do-nothing condition so that the pressure takes a specified value on the outlet - in a mean sense. The boundary conditions are of the form

\mu \partial_n u - pn = g

and g will be the mean value of the pressure on the boundary (up to a sign).

Please consider the book / paper I have linked for more details. I hope that this helps you solve your problem.

1 Like

Thank you for your reply, but I can use COMSOL to solve this problem under the Dirichlet BC on the pressure. Might be there are some details I don’t know.

Thank you for your help, and I will try your suggestion.

Hello, I tried your suggestion.

For my case, the normal vector on inlet is [-1, 0]^T and outlet is [1, 0]^T. If I’m right, the inlet condition and outlet conditon should be reformulated as follow:

\begin{cases}-p+2\mu\frac{\partial u}{\partial x} = -p_0,\text{ on }\Gamma_{inlet}\\-p+2\mu\frac{\partial u}{\partial x} = 0,\text{ on }\Gamma_{outlet}\\u_y=0\text{ on }\Gamma_{inlet}\cup\Gamma_{outlet}\end{cases}

Thus, I add the weak formulation term \int (p-\mu\frac{\partial u}{\partial x})q.

But the result is still wrong, the velocity field seems to be zero.

Here is my code

import dolfinx.fem.petsc
import gmsh
import numpy as np
import ufl
import mpi4py
from basix.ufl import element, mixed_element
from dolfinx import fem
from dolfinx.fem import Function, dirichletbc, form, functionspace, locate_dofs_topological
from ufl import div, dx, grad, inner, ds, dot, nabla_grad, sym, Identity, FacetNormal, Dx
import viskex
import dolfinx
from petsc4py import PETSc

# geometry & mesh parameter (unit is m)
L, meshsz = 0.001, 2E-5
WIDTH, HEIGHT = 5 * L, L

# fluid parameter
rho, mu = 1E3, 1E-3  # kg/m^3, Pa * s
init_pressure = 2  # Pa

gmsh.initialize()
gmsh.option.setNumber("Mesh.MshFileVersion", 2.2)
gmsh.option.setNumber("Mesh.Algorithm", 5)
gmsh.option.setNumber('General.Terminal', 0)
gmsh.option.setNumber('Mesh.MeshSizeMax', meshsz)

geo = gmsh.model.occ.addRectangle(-WIDTH/2, 0, 0, WIDTH*2, HEIGHT)
gmsh.model.occ.synchronize()
gmsh.model.mesh.generate(2)

bd = gmsh.model.occ.getEntities(1)
wall, inlet, outlet = [], [], []
for reg in bd:
    tag = reg[1]
    x, y, _ = gmsh.model.occ.getCenterOfMass(1, tag)
    if x < -WIDTH / 2 + meshsz / 2:
        inlet.append(tag)
    elif x > 3 * WIDTH / 2 - meshsz / 2:
        outlet.append(tag)
    else:
        wall.append(tag)
gmsh.model.addPhysicalGroup(1, inlet, 1, name="inlet")
gmsh.model.addPhysicalGroup(1, outlet, 2, name="outlet")
gmsh.model.addPhysicalGroup(1, wall, 3, name="wall")
gmsh.model.addPhysicalGroup(2, [geo], 0, name="solid")

mesh, subdomains, boundaries = dolfinx.io.gmshio.model_to_mesh(gmsh.model, comm=mpi4py.MPI.COMM_WORLD, rank=0, gdim=2)
mesh.topology.create_connectivity(mesh.topology.dim - 1, mesh.topology.dim)
inlet = boundaries.indices[boundaries.values == 1]
outlet = boundaries.indices[boundaries.values == 2]
wall = boundaries.indices[boundaries.values == 3]
viskex.dolfinx.plot_mesh_tags(mesh, boundaries, "boundaries")


# refer https://docs.fenicsproject.org/dolfinx/main/python/demos/demo_stokes.html
def mixed_direct(msh, inlet, outlet, wall):

    def inlet_pressure_expression(x):
        return np.ones(x.shape[1]) * init_pressure

    def inlet_velocity_expression(x):
        return np.stack((1 * np.ones(x.shape[1]), np.zeros(x.shape[1])))

    P2 = element("Lagrange", msh.basix_cell(), 2, shape=(msh.geometry.dim,))
    P1 = element("Lagrange", msh.basix_cell(), 1)

    TH = mixed_element([P2, P1])
    W = functionspace(msh, TH)
    W0, W1 = W.sub(0), W.sub(1)
    V, _ = W0.collapse()
    Q, _ = W1.collapse()

    # No slip boundary condition
    noslip = Function(V)
    dofs = locate_dofs_topological((W0, V), 1, wall)
    bc0 = dirichletbc(noslip, dofs, W0)

    # u \cdot t = 0 on inlet and outlet, t=[0, 1]^T
    W00 = W0.sub(1)
    V0, _ = W00.collapse()
    symmetry = Function(V0)
    dofs = locate_dofs_topological((W00, V0), 1, np.hstack([inlet, outlet]))
    bc1 = dirichletbc(symmetry, dofs, W00)

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

    # Define variational problem
    (u, p) = ufl.TrialFunctions(W)
    (v, q) = ufl.TestFunctions(W)
    f = Function(V)

    a = form(
        (mu * inner(grad(u), grad(v)) - inner(p, div(v)) + inner(div(u), q)) * dx
        + (inner(p, q) - 2 * mu * inner(u[0].dx(0), q)) * ds(1)
        + (inner(p, q) - 2 * mu * inner(u[0].dx(0), q)) * ds(2)
    )
    L = form(inner(f, v) * dx + init_pressure * q * ds(1))

    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")

    # Configure MUMPS to handle pressure nullspace
    pc = ksp.getPC()
    pc.setType("lu")
    pc.setFactorSolverType("mumps")
    pc.setFactorSetUpSolverType()
    pc.getFactorMatrix().setMumpsIcntl(icntl=24, ival=1)
    pc.getFactorMatrix().setMumpsIcntl(icntl=25, ival=0)

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

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

    viskex.dolfinx.plot_vector_field(u, "u", glyph_factor=1e-4)
    viskex.dolfinx.plot_scalar_field(p, "p")


mixed_direct(mesh, inlet=inlet, outlet=outlet, wall=wall)

You need to correctly impose the boundary conditions… The do-nothing boundary condition is a natural one, so that you have to impose these weakly. I made a MWE with FEniCS (not FEniCSx, I don’t have any experience with this).

from fenics import *
import matplotlib.pyplot as plt

mesh = UnitSquareMesh(8, 8)
dx = Measure("dx", domain=mesh)
n = FacetNormal(mesh)

boundaries = MeshFunction("size_t", mesh, 1)
left = CompiledSubDomain("on_boundary && near(x[0], 0.0)")
right = CompiledSubDomain("on_boundary && near(x[0], 1.0)")
no_slip = CompiledSubDomain("on_boundary && (near(x[1], 0.0) || near(x[1], 1.0))")
left.mark(boundaries, 1)
right.mark(boundaries, 2)
no_slip.mark(boundaries, 3)
ds = Measure("ds", domain=mesh, subdomain_data=boundaries)

u_elem = VectorElement("CG", mesh.ufl_cell(), 2)
p_elem = FiniteElement("CG", mesh.ufl_cell(), 1)
V = FunctionSpace(mesh, MixedElement(u_elem, p_elem))

u, p = TrialFunctions(V)
v, q = TestFunctions(V)

p_left = 2.0
p_right = -3.0
a = inner(grad(u), grad(v)) * dx - p * div(v) * dx - q * div(u) * dx
L = -Constant(p_left) * dot(n, v) * ds(1) - Constant(p_right) * dot(n, v) * ds(2)

bc_no_slip = DirichletBC(V.sub(0), Constant((0.0, 0.0)), no_slip)
bc_normal_left = DirichletBC(V.sub(0).sub(1), Constant(0.0), left)
bc_normal_right = DirichletBC(V.sub(0).sub(1), Constant(0.0), right)

bcs = [bc_no_slip, bc_normal_left, bc_normal_right]

up = Function(V)
solve(a == L, up, bcs)
u, p = up.split(True)

mean_pressure_left = assemble(p * ds(1)) / assemble(Constant(1.0) * ds(1))
mean_pressure_right = assemble(p * ds(2)) / assemble(Constant(1.0) * ds(2))

print(f"Prescribed p left: {p_left:.6e}  Mean p left: {mean_pressure_left:.6e}")
print(f"Prescribed p right: {p_right:.6e}  Mean p right: {mean_pressure_right:.6e}")

plt.figure()
plot(u)
plt.figure()
plot(p)
plt.show()

and it works:

Prescribed p left: 2.000000e+00  Mean p left: 2.000000e+00
Prescribed p right: -3.000000e+00  Mean p right: -3.000000e+00

Below are the velocity and pressure fields

discourse_u
discourse_p

Hope this helps

It seems you want to apply \boldsymbol{\sigma} \cdot \mathbf{n} = \mathbf{t}_{inlet} on \Gamma_{inlet} where \mathbf{t}_{inlet} is the given traction.

I have modified your code as follows which works on v0.9.0.
Hope you find it useful.

import dolfinx.fem.petsc
import numpy as np
import ufl
import mpi4py
from mpi4py import MPI
from basix.ufl import element, mixed_element
from dolfinx import fem
from dolfinx.mesh import (CellType, create_rectangle,
                          locate_entities_boundary, locate_entities, meshtags)
from dolfinx.fem import Function, dirichletbc, form, functionspace, locate_dofs_topological, Constant
from ufl import div, dx, grad, inner, ds, dot, nabla_grad, sym, Identity, FacetNormal, Dx, Measure
import viskex
import dolfinx
from petsc4py import PETSc

# fluid parameter
rho, mu = 1E3, 1E-3  # kg/m^3, Pa * s
traction = 2.0  # Pa

# Create mesh

msh = create_rectangle(MPI.COMM_WORLD,
                       [np.array([0, 0]), np.array([1, 1])],
                       [16, 16],
                       CellType.triangle)

# Function to mark walls

def boundary_noslip(x):
    return np.logical_or(np.isclose(x[1], 0.0),
                         np.isclose(x[1], 1.0))

# boundaries

boundaries = [(1, lambda x: np.isclose(x[0], 0)),
              (2, lambda x: np.isclose(x[0], 1)),
              (3, lambda x: np.isclose(x[1], 0)),
              (4, lambda x: np.isclose(x[1], 1))]

facet_indices, facet_markers = [], []
fdim = msh.topology.dim - 1
for (marker, locator) in boundaries:
    facets = locate_entities(msh, fdim, locator)
    facet_indices.append(facets)
    facet_markers.append(np.full_like(facets, marker))
facet_indices = np.hstack(facet_indices).astype(np.int32)
facet_markers = np.hstack(facet_markers).astype(np.int32)
sorted_facets = np.argsort(facet_indices)
facet_tag = meshtags(msh, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])

ds = Measure("ds", domain=msh, subdomain_data=facet_tag)

# function for noslip velocity

def noslip_expression(x):
    return np.stack((np.zeros(x.shape[1]), np.zeros(x.shape[1])))

def mixed_direct(msh):

    P2 = element("Lagrange", msh.basix_cell(), 2, shape=(msh.geometry.dim,))
    P1 = element("Lagrange", msh.basix_cell(), 1)

    TH = mixed_element([P2, P1])
    W = functionspace(msh, TH)
    W0, W1 = W.sub(0), W.sub(1)
    V, _ = W0.collapse()
    Q, _ = W1.collapse()

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

    # Collect Dirichlet boundary conditions
    bcs = [bc0]

    # Define variational problem
    (u, p) = ufl.TrialFunctions(W)
    (v, q) = ufl.TestFunctions(W)
    f = Constant(msh, (PETSc.ScalarType(0), PETSc.ScalarType(0)))  # type: ignore
    g = Constant(msh, (PETSc.ScalarType(traction), PETSc.ScalarType(0)))  # type: ignore

    a = form(
            (mu * inner(grad(u), grad(v)) - inner(p, div(v)) + inner(div(u), q)) * dx
            )

    L = form(
            inner(f, v) * dx + inner(g, v) * ds(1)
            )

    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")

    # Configure MUMPS to handle pressure nullspace
    pc = ksp.getPC()
    pc.setType("lu")
    pc.setFactorSolverType("mumps")
    pc.setFactorSetUpSolverType()
    pc.getFactorMatrix().setMumpsIcntl(icntl=24, ival=1)
    pc.getFactorMatrix().setMumpsIcntl(icntl=25, ival=0)

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

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

    viskex.dolfinx.plot_vector_field(u, "u", glyph_factor=1e-4)
    viskex.dolfinx.plot_scalar_field(p, "p")


mixed_direct(msh)

I want to add a comment that in this problem the pressure field is uniquely determined, so the following lines are not really necessary.

    pc.getFactorMatrix().setMumpsIcntl(icntl=24, ival=1)
    pc.getFactorMatrix().setMumpsIcntl(icntl=25, ival=0)