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.
