Stokes Operator is Singular After Applying Dirichlet 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)
2

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”

Answered here: `addboundaryfacetset!` adding inconsistant number of facets · Issue #1273 · Ferrite-FEM/Ferrite.jl · GitHub, but for others seeing this:
probably need to use all = false in addboundaryfacetset! for "boundary_other"

3

Not sure.

Is inadvertedly a boundary condition on the pressure field missing, causing the pressure block to be singular?

4

Even if I fix the pressure field, the operator remains singular. Since the pressure is defined up to addition by a constant, I have fixed the pressure field using

addnodeset!(
    grid, "p_anchor", 
    x -> x ≈ grid.nodes[1].x
)

add!(ch_ns, Dirichlet(
    :p, getnodeset(grid, "p_anchor"), x -> 1.0
))
5

Right.

Again, not sure.

Is this error caused by nuisances in the mesh?

Should one try a 2D rectangular or 3D mesh cube mesh using the built-in Ferrite mesh generator? See e.g. finite_element_electrical_engineering/project-based-assignment/metal-hydride-storage/notebooks/stokes_channel.ipynb at main · ziolai/finite_element_electrical_engineering · GitHub

Good luck.

6

Yes, you are correct. Something is wrong with my mesh. I created a cylinder with a sphere removed from its center as below, adjusted cyl_wall_pred, and now everything is working as expected. Any suggestions for how I can diagnose/fix the mesh? I have the original model in FreeCAD, and used gmsh to create the .unv file.

Here is how I generated the test grid:

function create_grid()
    gmsh.initialize()
    gmsh.clear()

    R_cyl = 1.0
    H_cyl = 5.0

    R_sph = 0.5

    lc = 0.2   # target mesh size

    cyl = gmsh.model.occ.addCylinder(
        0.0, 0.0, 0.0,     # base center
        0.0, 0.0, H_cyl,  # axis direction
        R_cyl
    )

    sph = gmsh.model.occ.addSphere(
        0.0, 0.0, H_cyl/2,
        R_sph
    )

    cut_entities, _ = gmsh.model.occ.cut(
        [(3, cyl)],
        [(3, sph)]
    )

    gmsh.model.occ.synchronize()

    gmsh.model.mesh.generate(3)
    
    grid = togrid()

    gmsh.finalize()

    return grid
end

grid = create_grid()
7

Nice!

Is it beneficial here to elaborate on how you define the boundary patches?

Is it wise here to perform tests in 2D first?

8

I have recreated the mesh using a different Stanford bunny model and everything is now working as expected! I do not know what the issue was with the original mesh. For anyone interested, you can find the new mesh here.

Thanks for your help!

1 Like