Hi,
I am pretty new to all this - FEM and dolfinx.
I am trying to get the frequency response of an acoustic cavity using modal superposition. Essentially, I am follwoing this tutorial: Modal Superposition in Acoustic Cavities – Unda Project
Although my knowledge is pretty limited, I was able to follow the tutorials on that site, with lots of trying around and referring to other resources.
But I really can’t figure out the last step that in this tutorial is described as “a couple of matrix multiplications” and computing the modal participation factors using a ksp linear solver.
I’ve searched a lot online and even gave chatgpt a shot. The fact that dolfinx and dolfin seem to have quite some differences doesnt make it easier for a beginner to try out code found online.
Any help is greatly appreciated.
Below is the code i have currently.
It seems to correctly extract the modes of a half-open tube. The open-ness is modeled with a dirichlet bc at the right.
I would be interested to apply an acoustic load (As far as I read it should be a velocity source, but i don’t quite know) to the left of the tube and see what the frequency response looks like via modal superposition.
Thank you!
# Dolfinx version: 0.6.0, complex number support.
# Python Version: 3.10.9
#
# IMPORTS
import matplotlib.pyplot as plt
from dolfinx.mesh import create_box
import pyvista as pv
pv.set_jupyter_backend('static')
print("Pyvista Backend: ", pv.global_theme.jupyter_backend)
from dolfinx import plot
from petsc4py import PETSc
import numpy as np
import dolfinx
import numpy as np
import ufl
import sys, slepc4py
slepc4py.init(sys.argv)
from dolfinx import fem
from dolfinx.fem import (Constant, dirichletbc, Function, FunctionSpace, VectorFunctionSpace,
locate_dofs_topological)
from dolfinx.mesh import (create_box, meshtags, locate_entities, CellType,
locate_entities_boundary)
from dolfinx.io import XDMFFile
from mpi4py import MPI
from petsc4py import PETSc
from slepc4py import SLEPc
from dolfinx.mesh import create_rectangle
from ufl import (inner, grad, dx)
from dolfinx.fem import Function, FunctionSpace, assemble_scalar, form, Constant
from tqdm import tqdm
# -----------------------------------------
deg = 1
c0 = 340 # [m/s]
freqMax = 20000 #Hz Maximum Frequency of interest
λ_min = c0/freqMax
meshSizeMin = λ_min/10. # [m] rule of thumb according to https://undaproject.com/sound-propagation-in-3d-cavities
print(f" - Speed of sound: {c0} m/s")
print(f' - Lambda min: {λ_min} m')
print(f" - minimum Mesh Size approx: {meshSizeMin*1e3:.2f} mm")
L_mm = 500 # [mm] Pipe Length
D_mm = 20 # [mm] Pipe Diameter
L = L_mm/1e3
D = D_mm/1e3
f_analytic = c0/(4*L)
print(f"Predicted Fundamental Frequency: {f_analytic} Hz")
print("Expected Harmonics:")
for i in range(5):
print(f"{i+1}. { f_analytic*((i+1)*2+1)} Hz")
print('....')
## Deinition of Geometry
ll = np.array([0,0]) #lower left corner [mm]
ur = np.array([L_mm,D_mm]) #upper right corner [mm]
ll = ll/1e3 #to meters
ur = ur/1e3 #to meters
dims = ur - ll
numCellsPerDim = np.ceil(dims / meshSizeMin).astype(np.int64)
print(f"N Cells: {numCellsPerDim}")
estCellsTotal = int(numCellsPerDim[0]*numCellsPerDim[1])
print(f'Estimated Cells Total: {estCellsTotal}')
## Mesh Creation
mesh = create_rectangle(MPI.COMM_WORLD, [list(ll), list(ur)], list(numCellsPerDim), cell_type=dolfinx.mesh.CellType.quadrilateral)
tdim = mesh.topology.dim
#pv.start_xvfb()
topology, cell_types, geom = plot.create_vtk_mesh(mesh, tdim)
grid = pv.UnstructuredGrid(topology, cell_types, geom)
plotter = pv.Plotter(window_size=[1500, 500])
plotter.add_mesh(grid, show_edges=True)
plotter.view_xy()
plotter.show_grid(font_size = 15)
plotter.camera.zoom(2.5)
plotter.show()
# Test and trial function space
V = FunctionSpace(mesh, ("CG", deg))
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
# Define boundary conditions
def getRightBorder(x):
return np.isclose(x[0], ur[0])
fdim = mesh.topology.dim - 1
# -----------------------------------------------------------------------
boundary_facetsR = locate_entities_boundary(mesh, fdim, getRightBorder)
print("Boundary facets right:", boundary_facetsR)
u_D = Function(V)
u_D.interpolate(lambda x: x[1]*0) #Doesnt make much sense to me but works
bcR = dirichletbc(u_D, locate_dofs_topological(V, fdim, boundary_facetsR))
# Stiffness and Mass Matrices
stiff_matrix = inner(grad(u), grad(v)) * dx
mass_matrix = inner(u, v) * dx
# Using the "diagonal" kwarg ensures that Dirichlet BC modes will not be among
# the lowest-frequency modes.
SM = dolfinx.fem.petsc.assemble_matrix(form(stiff_matrix), [bcR], diagonal=62831)
MM = dolfinx.fem.petsc.assemble_matrix(form(mass_matrix), [bcR], diagonal=1/62831)
SM.assemble()
MM.assemble()
# Setup the eigensolver
N = 50 #num modes to extract
problem = SLEPc.EPS().create()
A = [ ]
A.append(SM)
A.append(MM)
# Setup Spectral Transformation
st = SLEPc.ST().create()
st.setType(SLEPc.ST.Type.SINVERT)
st.setShift(0.1)
st.setFromOptions()
problem.setST(st)
problem.setDimensions(N)
problem.setProblemType(SLEPc.EPS.ProblemType.GHEP)
problem.setInterval(10,18)
problem.setOperators(SM,MM)
# Solve the eigenproblem
problem.solve()
# Getting the result vectors
xr, xi = SM.createVecs()
nconv = problem.getConverged()
if nconv > 0:
modes = []
eig_values = []
eig_vectors = []
eig_freq = []
vrs = []
for i in range(nconv):
# eigenvector = Function(V)
# modes.append(eigenvector)
k = problem.getEigenpair(i, xr, xi)
fn = np.sqrt(k.real)/(2*np.pi)*c0
eig_values.append(k)
eig_freq.append(fn)
vrs.append(xr.array.copy())