Matrix structure of Stokes evp FeniCSx

dolfinx
1

Hello. I’m trying to solve the stokes generalized eigenvalue problem using Taylor-Hood elements. When i try to spy the matrices to see their structure, i get this diagonal pattern

A_csr
A_csr434×418 8.56 KB

(M is similar to this one)
I was expecting a different pattern at least in the theory behind it, does this have something to do with how fenicsx handles the dofs or the manner it solves the eigenvalue problem?. I’m new to dolfinx and FEM in general.

Here’s my code, thanks

from mpi4py import MPI
import numpy as np
from scipy.sparse import csr_matrix
import matplotlib.pylab as plt
import dolfinx 

import ufl
from basix.ufl import element, mixed_element
from dolfinx import default_real_type
from dolfinx.fem import (
    Function,
    dirichletbc,
    form,
    functionspace,
    locate_dofs_topological,
)
from dolfinx.fem.petsc import assemble_matrix
from dolfinx.mesh import CellType, create_rectangle, locate_entities_boundary
from ufl import div, dx, grad, inner


# Mesh & elements
msh = create_rectangle(MPI.COMM_WORLD, [np.array([0,0]), np.array([1,1])],  [16,16],  CellType.triangle)
P2 = element("Lagrange", msh.basix_cell(), 2, shape=(msh.geometry.dim,), dtype=default_real_type)
P1 = element("Lagrange", msh.basix_cell(), 1, dtype=default_real_type)
TH = mixed_element([P2, P1])
W = functionspace(msh, TH)

def noslip_boundary(x):
    return np.isclose(x[0], 0.0) | np.isclose(x[0], 1.0) | np.isclose(x[1], 0.0) | np.isclose(x[1], 1.0)


# no-slip boundary condition
W0 = W.sub(0) 
V0, _ = W0.collapse() 

noslip = Function(V0) 
facets = locate_entities_boundary(msh, 1, noslip_boundary)
dofs = locate_dofs_topological((W0, V0), 1, facets)
bc = dirichletbc(noslip, dofs, W0)


# Variational formulation
(u,p) = ufl.TrialFunctions(W) 
(v,q) = ufl.TestFunctions(W)

a = form((inner(grad(u), grad(v)) - inner(div(v), p) + inner(div(u), q))* dx )
m = form(inner(u,v)* dx)

A = assemble_matrix(a, bcs=[bc])
A.assemble()

M = assemble_matrix(m, bcs=[bc])
M.assemble()


# Get CSR representation
I, J, V = A.getValuesCSR()
A_csr = csr_matrix((V, J, I))
plt.spy(A_csr)
plt.show()
2

Which pattern were you expecting?

3

My error, i may have attached it
This is the pattern i expeted. This is from a friend that implemented it in firedrake

imagen
imagen968×510 79 KB

4

Perhaps you’re expecting a block matrix? If you assemble the system as a block matrix system (e.g., injecting your code at roughly this point in the demo), you should see something like:

image
image770×771 19 KB

Bear in mind that DOLFINx reorders DoFs to minimise matrix bandwidth which is important for computational efficiency.

1 Like
5

Thanks a lot!, that solves my question