DOLFINx/SLEPc eigenfrequency problem: spurious modes and frequency scaling issue

dolfinx
1

Hi folks,

I’m using the SLEPc eigensolver with DOLFINx for the first time, and I can’t seem to get correct results even for a simple benchmark problem. For context, I’m running the Python version from conda, installed via

conda create -n fenicsx-env -c conda-forge fenics-dolfinx matplotlib ipython

The problem is an eigenmode analysis of a beam of size 0.3 \times 0.05 \times 0.02 m, fixed at one end, with parameters E = 1\times10^6, \nu = 0.3, and \rho = 1000. A simple analytical calculation, as well as a COMSOL simulation, gives the first few natural frequencies around 1.2 Hz, 2.8 Hz, etc.

However, when I run the MWE attached below, the code produces the first two modes with eigenvalues exactly equal to 1. Starting from the third mode, the results begin to look reasonable, but the eigenvalues are still much too large. In fact, the FEniCS results appear to follow the correct sequence of modes, but all frequencies are about 6.3× larger than expected (e.g. 7.5 Hz, 17.7 Hz, …). The first two modes (with eigenvalue 1) also look strange: either oscillations localized near the Dirichlet boundary or some sort of whole-body twisting. I can’t upload images but the MWE can generate those modes.

This makes me suspect I don’t have the problem or the boundary conditions set up correctly, but I can’t see where the mistake is. Any advice or insight would be greatly appreciated!

MWE:

from mpi4py import MPI
from petsc4py import PETSc
from slepc4py import SLEPc
from petsc4py.PETSc import ScalarType

import numpy as np
import ufl
from basix.ufl import element
import os

from dolfinx import fem, io, mesh
from dolfinx.mesh import locate_entities_boundary, exterior_facet_indices
from dolfinx.fem.petsc import assemble_matrix

# 1. Create 3D mesh and function space
# Domain dimensions
Lx, Ly, Lz = 0.3, 0.05, 0.02   # box size
nx, ny, nz = 30, 15, 5      # number of elements in each direction
# Create a structured box mesh
domain = mesh.create_box(
    MPI.COMM_WORLD,
    [[0.0, 0.0, 0.0], [Lx, Ly, Lz]],   # two diagonal corners
    [nx, ny, nz],
    cell_type=mesh.CellType.hexahedron
)

# 2. Material parameters
E, nu, rho = 1e6, 0.3, 1000.
mu = E / (2 * (1 + nu))
lmbda = E * nu / ((1 + nu) * (1 - 2 * nu))

P1 = element("Lagrange", domain.basix_cell(), degree=1, shape=(3,))
V = fem.functionspace(domain, P1)
print(f"Value size of function space: {V.dofmap.index_map_bs}")

# 3. Define strain and stress
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)

def epsilon(u): return ufl.sym(ufl.grad(u))
def sigma(u): return lmbda * ufl.tr(epsilon(u)) * ufl.Identity(3) + 2 * mu * epsilon(u)

# 4. Bilinear forms
a_form = ufl.inner(sigma(u), epsilon(v)) * ufl.dx
m_form = rho * ufl.inner(u, v) * ufl.dx

# 5. Apply Dirichlet BC (fix boundary faces where all facet nodes have x^2 + z^2 < 25)
# Identify exterior facets (codim=1 = dim-1 = 2 for 3D)
boundary_facets = locate_entities_boundary(domain, domain.topology.dim - 1,
    lambda x: np.isclose(x[0], 0.0))

# Apply BC on all components
bcs = []
for i in range(3):  # x, y, z components
    dofs = fem.locate_dofs_topological(V.sub(i), domain.topology.dim - 1, boundary_facets)
    bc = fem.dirichletbc(ScalarType(0), dofs, V.sub(i))
    bcs.append(bc)

# 6. Assemble stiffness and mass matrices
a = fem.form(a_form)
m = fem.form(m_form)
A = assemble_matrix(a, bcs=bcs)
M = assemble_matrix(m, bcs=bcs)
A.assemble()
M.assemble()

# print(f"||A||_F = {A.norm()}, ||M||_F = {M.norm()}")

# 7. Solve eigenvalue problem
eps = SLEPc.EPS().create(MPI.COMM_WORLD)
eps.setOperators(A, M)
eps.setProblemType(SLEPc.EPS.ProblemType.GHEP)
eps.setType(SLEPc.EPS.Type.KRYLOVSCHUR)
eps.setDimensions(nev=50, ncv=100)
eps.setTolerances(tol=1e-4, max_it=10000)
eps.setWhichEigenpairs(SLEPc.EPS.Which.SMALLEST_REAL)
st = eps.getST()
ksp = st.getKSP()
pc = ksp.getPC()
pc.setType("ilu")     

eps.solve()

n_conv = eps.getConverged()
if domain.comm.rank == 0:
    print(f"Number of converged eigenpairs: {n_conv}")
    for i in range(min(n_conv, 10)):
        eigval = eps.getEigenvalue(i)
        omega = np.sqrt(eigval.real) if eigval.real > 0 else 0.0
        print(f"Mode {i}: λ = {eigval:.6e}, ω = sqrt(λ) = {omega:.6e}")

# 8. Export first few modes to file
u_mode = fem.Function(V)
vec = PETSc.Vec().create(MPI.COMM_WORLD)
vec.setSizes(V.dofmap.index_map.size_global * V.dofmap.index_map_bs)
vec.setFromOptions()
vec.assemble()

needed_modes = min(n_conv, 50)
with io.VTKFile(domain.comm, "eigenmodes_3D.pvd", "w") as vtk:
    for i in range(needed_modes):
        eigval = eps.getEigenvalue(i)
        eps.getEigenpair(i, vec)
        u_mode.x.array[:] = vec.array
        vtk.write_function(u_mode, i)
        np.save(f"mode_{i}.npy", u_mode.x.array)
        with open("eigenvals.txt", "a") as f:
            f.write(f"{i} {eps.getEigenvalue(i).real}\n")
2

Be careful, the FEniCSx modes don’t result 6.3 times larger, but exactly 2\pi times larger. You are confusing the frequency f [Hz] and angular frequency \omega [rad/s], linked by the relation:

\begin{equation} \omega = 2\pi f \end{equation}

The first two are spurious modes, caused by the lifting operation used to enforce the Dirichlet boundary condition.
One way to remove them is to assign a very large stiffness to the corresponding boundary degrees of freedom.
Another option is simply to ignore them, since they are non-physical.

1 Like
3

I see. This resolves the problem and thank you!