Solving eigenmodes of a 2D rectangular cavity

Greetings everyone,
I am fairly new to Fenicsx, I did the tutorials available in the Fenicsx documentation, read some examples but that’s about it.

I want to study the acoustic eigenmodes of an acoustic cavity with Fenicsx, to start easily, I tried doing it on a simple 2D rectangular cavity. I tried applying the same method as in the tutorial by setting up Dirichlet boundary condtions and using the Helmholtz equation. However, I am facing a problem when it comes to defining the equation. Since I want to study the eigen modes of the cavity, I don’t know the value of the wave number in the helmholtz equations. For the moment I tried implementing by adapting this :
https://docs.fenicsproject.org/dolfinx/v0.7.2/python/demos/demo_helmholtz.html

but without much success :
mport numpy as np

import ufl
from dolfinx.fem import Function, FunctionSpace
from dolfinx.fem.petsc import LinearProblem
from dolfinx.io import XDMFFile
from ufl import dx, grad, inner
import dolfinx
from mpi4py import MPI

wavenumber still need to be given

k0 = 2.75 * np.pi

L=1

number of elements in each direction of msh

n_elem = 128

msh = mesh = dolfinx.mesh.create_rectangle(
MPI.COMM_WORLD,
[[0, 0], [L, L]],
[n_elem, n_elem])
n = ufl.FacetNormal(msh)

Source amplitude

A = 1

Test and trial function space

V = FunctionSpace(msh, (“Lagrange”, 1))

Define variational problem

u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
f = Function(V)
f.interpolate(lambda x: A * k02 * np.cos(k0 * x[0]) * np.cos(k0 * x[1]))
a = inner(grad(u), grad(v)) * dx - k0
2 * inner(u, v) * dx
L = inner(f, v) * dx

#Implementing boudary conditons

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

def BXL(x):
return np.isclose(x[0],L)

def BY0(x):
return np.isclose(x[1],0)

def BYL(x):
return np.isclose(x[1],L)

bx0_dofs = dolfinx.fem.locate_dofs_geometrical(V, BX0)
bxl_dofs = dolfinx.fem.locate_dofs_geometrical(V, BXL)
by0_dofs = dolfinx.fem.locate_dofs_geometrical(V, BY0)
byl_dofs = dolfinx.fem.locate_dofs_geometrical(V, BYL)

bcxo= dolfinx.fem.dirichlet(dolfinx.default_scalar_type(0), bx0_dofs, V)
bcxl= dolfinx.fem.dirichlet(dolfinx.default_scalar_type(0), bxl_dofs, V)
bcyo= dolfinx.fem.dirichlet(dolfinx.default_scalar_type(0), by0_dofs, V)
bcyl= dolfinx.fem.dirichlet(dolfinx.default_scalar_type(0), byl_dofs, V)

bcs = [bcxo, bcxl, bcyo, bcyl]

Compute solution

uh = Function(V)
uh.name = “u”
problem = LinearProblem(a, L, u=uh,bcs=bcs, petsc_options={“ksp_type”: “preonly”, “pc_type”: “lu”})
problem.solve()

Save solution in XDMF format (to be viewed in Paraview, for example)

with XDMFFile(MPI.COMM_WORLD, “/home/daep/mar.baudry/Documents/results/plane_wave.xdmf”, “w”, encoding=XDMFFile.Encoding.HDF5) as file:
file.write_mesh(msh)
file.write_function(uh)

Is there a way to define a PDE with an unknown parameter or maybe is there another method ? I also tried using Fenics eigenvalues solver but I am not convinced of the result and I’m not sure I understand what the matrix are representing :

from dolfin import *
import numpy as np
import fenics
import ufl

Test for PETSc and SLEPc

if not has_linear_algebra_backend(“PETSc”):
print (“DOLFIN has not been configured with PETSc. Exiting.”)
exit()

if not has_slepc():
print (“DOLFIN has not been configured with SLEPc. Exiting.”)

exit()

Define mesh, function space

lx = 8
ly = 8
c = 341
f = 57.17 # Frequency of the Simulation, Hz
s = c / (10 * f) # Element Size
Nx = np.intp(np.ceil(lx / s)) # Number of elements for each direction
Ny = np.intp(np.ceil(ly / s))
tol = 1e-10 # Tolerance for boundary condition definitions

mesh = UnitSquareMesh(lx, ly)

V = FunctionSpace(mesh, “Lagrange”, 1)

Define basis and bilinear form

u = TrialFunction(V)
v = TestFunction(V)

#a = inner(grad(u), grad(v)) * dx - k_sq * inner(u, v) * dx
#L = inner(f, v) * dx
a = (ufl.inner(ufl.nabla_grad(u), ufl.nabla_grad(v)) * ufl.dx)
m = inner(u, v) * dx

Assemble stiffness form

A = PETScMatrix()
assemble(a, tensor=A)

M = PETScMatrix()
assemble(m, tensor=M)

#Boundary conditions
def BX0(x, on_boundary):
return on_boundary and fenics.near(x[0], 0, tol)

def BXL(x, on_boundary):
return on_boundary and fenics.near(x[0], lx, tol)

def BY0(x, on_boundary):
return on_boundary and fenics.near(x[1], 0, tol)

def BYL(x, on_boundary):
return on_boundary and fenics.near(x[1], ly, tol)

#Definition of boundary conditions here for a completely closed box

bcx0 = DirichletBC(V, 0, BX0)
bcxl = DirichletBC(V, 0, BXL)
bcy0 = DirichletBC(V, 0, BY0)
bcyl = DirichletBC(V, 0, BYL)

bcs = [bcx0, bcxl, bcy0, bcyl]

bcx0.apply(A)
bcx0.apply(M)

bcxl.apply(A)
bcxl.apply(M)

bcy0.apply(A)
bcy0.apply(M)

bcy0.apply(A)
bcyl.apply(M)

Create eigensolver

eigensolver = SLEPcEigenSolver(A,M)

Compute all eigenvalues of A x = \lambda x

print (“Computing eigenvalues. This can take a minute.”)
eigensolver.solve()

Extract largest (first) eigenpair

r, c, rx, cx = eigensolver.get_eigenpair(0)

print ("Largest eigenvalue: "), r

Initialize function and assign eigenvector

u = Function(V)
u.vector()[:] = rx

Plot eigenfunction

plot(u)
interactive()

Thanks for reading this and for your help !

Have a nice day

If you are willing to use non legacy dolfin but dolfinx instead, I’ve made a video on how to perform acoustic modal analysis:

https://www.youtube.com/watch?v=SH4yjPiGxXA

Best,

Antonio

2 Likes

Thanks a lot a will watch it with great interest !