Stokes Operator is Singular After Applying Drichlet Boundary Conditions

Specific Domains Modelling & Simulations
, ,
1

Hello! I am trying to implement a simple incompressible Navier-Stokes simulation. I have a mesh consisting of a cylinder with the Stanford bunny excised from its center. When I apply Drichlet boundary conditions, the Stokes operator becomes singular, and I cannot figure out why. I did not include it in the below example, but the problem persists even if I fix one of the pressure DOFs. You can find the mesh I am using here.

In case it helps you help me, I am well versed in linear algebra but not in how it applies to finite element analysis. In particular I am not familiar with all the lingo.

I have attached a minimal Jupyter notebook demonstrating my issue as well as the mesh.

using Ferrite
using LinearAlgebra
using SparseArrays, BlockArrays
using Ferrite
using FerriteGmsh
using FerriteGmsh: Gmsh

# load mesh
grid = togrid("/workspaces/BasisSubset/examples/advection-diffusion/bunny-cylinder.unv")

# specify velocity shape functions
ipu = Lagrange{RefTetrahedron, 2}()^3
qr = QuadratureRule{RefTetrahedron}(2)
cvu = CellValues(qr, ipu);

# specify pressure shape functions
ipp = Lagrange{RefTetrahedron, 1}() 
cvp = CellValues(qr, ipp);

# DOF handler for Navier-Stokes
dh_ns = DofHandler(grid)
add!(dh_ns, :u, ipu)
add!(dh_ns, :p, ipp)
close!(dh_ns);

# constraint handler for Navier-Stokes
ch_ns = ConstraintHandler(dh_ns);

# delete existing facesets
empty!(grid.facetsets)
empty!(grid.nodesets)

# determines if boundary is on wall of cylinder
cyl_wall_pred = x -> (abs(sqrt((x[1]+2.0)^2 + (x[2]+2.0)^2) - 50.0) < 1e-12 )

topo = ExclusiveTopology(grid)

# add facetset for everything on boundary
addboundaryfacetset!(
    grid, topo, "boundary_all",
    x -> true
)

# add facetset for everything on the wall of the cylinder
addboundaryfacetset!(
    grid, topo, "cylinder_wall", 
    x -> cyl_wall_pred(x)
)

# add facetset for everything that is not on the wall of the cylinder
addboundaryfacetset!(
    grid, topo, "boundary_other",
    x -> !cyl_wall_pred(x)
)

# WHAT I WANT
# add!(ch_ns, Dirichlet(
#     :u, getfacetset(grid, "cylinder_wall"), v -> begin
#         return (-v[2], v[1], 0)
#     end
# ))
# add!(ch_ns, Dirichlet(
#     :u, getfacetset(grid, "boundary_other"), v -> begin
#         return (0.0, 0.0, 0.0)
#     end
# ))

# BUT THIS DOES NOT WORK EITHER
add!(ch_ns, Dirichlet(
    :u, getfacetset(grid, "boundary_all"), v -> begin
        return (0.0, 0.0, 0.0)
    end
))

close!(ch_ns)
;

# maybe not relevant, but I noticed that the
# number of dofs in the facetsets "cylinder_wall"
# and "boundary_other" do not sum to that of 
# "boundary_all"
println(length(getfacetset(grid, "cylinder_wall")))
println(length(getfacetset(grid, "boundary_other")))
println(length(getfacetset(grid, "boundary_all")))

function assemble_stokes_operator!(K::SparseMatrixCSC, dh::DofHandler, cvu::CellValues, cvp::CellValues)
    n_basefuncs_u = getnbasefunctions(cvu)
    n_basefuncs_p = getnbasefunctions(cvp)
    n_basefuncs = n_basefuncs_u + n_basefuncs_p
    
    # will store element assembelies
    Ke = BlockedArray(
        Matrix{Float64}(undef, (n_basefuncs, n_basefuncs)), 
        [n_basefuncs_u, n_basefuncs_p], 
        [n_basefuncs_u, n_basefuncs_p]
    )

    # create assembler
    assembler_K = start_assemble(K)

    # loop over all cells
    for cell in CellIterator(dh)
        # update shape functions
        reinit!(cvu, cell)
        reinit!(cvp, cell)

        # clear element assembelies
        fill!(Ke, 0)

        # loop over quadrature points
        for q_point in 1:getnquadpoints(cvu)
            # get shape function coordinate transform
            dΩ = getdetJdV(cvu, q_point)

            # assemble velocity-velocity blocks
            for i in 1:n_basefuncs_u
                ∇a = shape_gradient(cvu, q_point, i)
                for j in 1:n_basefuncs_u
                    ∇u = shape_gradient(cvu, q_point, j)
                    Ke[BlockIndex((1,1), (i, j))] -= ∇a ⊡ ∇u * dΩ
                end
            end
            
            # assemble velocity-pressure blocks
            for j in 1:n_basefuncs_p
                p = shape_value(cvp, q_point, j)
                for i in 1:n_basefuncs_u
                    divu = shape_divergence(cvu, q_point, i)
                    Ke[BlockIndex((1,2), (i, j))] += (divu * p) * dΩ
                    Ke[BlockIndex((2,1), (j, i))] += (p * divu) * dΩ
                end
            end
        end

        # Assemble Ke into K
        assemble!(assembler_K, celldofs(cell), Ke)
    end
end

# allocate sparse matrices
K = allocate_matrix(dh_ns, ch_ns)

# allocate right hand side of solve
rhs = zeros(ndofs(dh_ns))

# asseble Stokes operator
assemble_stokes_operator!(K, dh_ns, cvu, cvp)

println("logdet before `apply!`: ", logdet(K))

# solve for velocity field
apply!(K, rhs, ch_ns)

println("logdet after `apply!`: ", logdet(K))

x0 = K \ rhs # this lines fails because K is singular
apply!(x0, ch_ns)