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 - k02 * 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