Mixed formulation of stokes equations

dolfinx
6

The issue is the boundary conditions.
The last argument of dirichletbc should not be V_sub, but W.sub(0) to indicate to the assemblers that it is the space we want to apply u_in on. I also changed the sign of your divergence term for symmetry.
MWE:

from mpi4py import MPI
from petsc4py import PETSc
import dolfinx as dfx
from dolfinx.fem.petsc import LinearProblem
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 noslip(x):
    return np.stack((np.zeros(x.shape[1]), np.zeros(x.shape[1])))

def inlet(x):
    return np.isclose(x[0], 0.0)

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

def outlet(x):
    return np.isclose(x[0], 5.0)

mesh = dfx.mesh.create_rectangle(MPI.COMM_WORLD, \
    [np.array([0.0, 0.0]), np.array([5, 1])], [50, 10], dfx.mesh.CellType.triangle)

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

ID_Wall = dfx.mesh.locate_entities_boundary(mesh, fdim, walls)
b_dofs_Wall = dfx.fem.locate_dofs_topological((W.sub(0),V_sub), fdim, ID_Wall)
u_wall = dfx.fem.Function(V_sub)
u_wall.interpolate(noslip)
bc_wall = dfx.fem.dirichletbc(u_wall, b_dofs_Wall, W.sub(0))

ID_Inlet = dfx.mesh.locate_entities_boundary(mesh, fdim, inlet)
b_dofs_Inlet = dfx.fem.locate_dofs_topological((W.sub(0),V_sub), fdim, ID_Inlet)
u_in = dfx.fem.Function(V_sub)
u_in.interpolate(inflow)
bc_inlet = dfx.fem.dirichletbc(u_in, b_dofs_Inlet, W.sub(0))

bc = [bc_inlet, bc_wall]

(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()
1 Like