I am working towards developing a monolithic stokes equation solver using mixed pressure velocity formulation and Taylor Hood elements. Having reviewed all the different posts, examples, and forums, I am able to get my implementation to run but get wrong results. Here is my attempt at replicating an MWE for flow through a long channel (aiming to replicate Poiseuille flow). At this point, I will highly benefit from someone to look at this and let me know why I am getting wrong results.
from mpi4py import MPI
from petsc4py import PETSc
import dolfinx as dfx
from dolfinx.fem.petsc import LinearProblem
from dolfinx.io import gmshio
import ufl as ufl
import numpy as np
def inflow(x):
return np.stack((np.full(x.shape[1], 1.0), np.zeros(x.shape[1])))
def wall(x):
return np.stack((np.zeros(x.shape[1]), np.zeros(x.shape[1])))
meshFile = 'test.msh'
ID_L = 5
ID_T = 6
ID_B = 8
mesh, c_id, f_id = gmshio.read_from_msh(meshFile, MPI.COMM_WORLD, gdim=2)
tdim = mesh.topology.dim
fdim = tdim - 1
V = ufl.VectorElement("Lagrange", mesh.ufl_cell(), 2)
P = ufl.FiniteElement("Lagrange", mesh.ufl_cell(), 1)
M = ufl.MixedElement([V,P])
W = dfx.fem.FunctionSpace(mesh, M)
V_sub, _ = W.sub(0).collapse()
P_sub, _ = W.sub(1).collapse()
b_dofs_L = dfx.fem.locate_dofs_topological((W.sub(0),V_sub), fdim, f_id.find(ID_L))
b_dofs_T = dfx.fem.locate_dofs_topological((W.sub(0),V_sub), fdim, f_id.find(ID_T))
b_dofs_B = dfx.fem.locate_dofs_topological((W.sub(0),V_sub), fdim, f_id.find(ID_B))
u_in = dfx.fem.Function(V_sub)
u_in.interpolate(inflow)
u_wall_1 = dfx.fem.Function(V_sub)
u_wall_1.interpolate(wall)
u_wall_2 = dfx.fem.Function(V_sub)
u_wall_2.interpolate(wall)
bc_1 = dfx.fem.dirichletbc(u_in, b_dofs_L, V_sub)
bc_2 = dfx.fem.dirichletbc(u_wall_1, b_dofs_T, V_sub)
bc_3 = dfx.fem.dirichletbc(u_wall_2, b_dofs_B, V_sub)
bc = [bc_1, bc_2, bc_3]
(us,ps) = ufl.TrialFunctions(W)
(vt,qt) = ufl.TestFunctions(W)
mu = dfx.fem.Constant(mesh, dfx.default_scalar_type(1.0))
b = dfx.fem.Constant(mesh, PETSc.ScalarType((0.0,0.0)))
weakForm = mu * ufl.inner(ufl.grad(vt), ufl.grad(us)) * ufl.dx \
- ps * ufl.div(vt) * ufl.dx \
+ qt * ufl.div(us) * ufl.dx \
- ufl.dot(vt,b) * ufl.dx
a = ufl.lhs(weakForm)
L = ufl.rhs(weakForm)
uh = dfx.fem.Function(W)
problem = LinearProblem(a, L, bcs=bc, u=uh)
problem.solve()
vtkFile = dfx.io.VTKFile(MPI.COMM_WORLD, "./mwe.pvd", "w")
vtkFile.write_function([uh.split()[0].collapse()])
vtkFile.close()
Happy to share more if needed.