Hello everyone,
I don’t know in which category I should put this, sorry. I just want to share this working code.
The following code solves the 2D Helmholtz equation for a square domain, places a monopole in the center of the square (you can move it, anyway) with flat volume velocity and computes the sound pressure - frequency spectrum in a specified location (microphone).
The boundaries are considered rigid wall (default BC).
In the end, saves a .csv file which columns are organized as follows:
[frequency real_part(pressure) imag_part(pressure]
and, finally plots the spectrum.
I am sorry if there are mistakes or “bad habits” in writing the code, I just started using FEniCSx 3 weeks ago, so I am still learning.
I am sure that people that work with Acoustic Simulation using commercial software like me will find this useful.
This is the scheme:
Code:
import numpy as np
import ufl
from dolfinx import geometry
from dolfinx.fem import Function, FunctionSpace, assemble_scalar, form, Constant
from dolfinx.fem.petsc import LinearProblem
from dolfinx.io import XDMFFile
from dolfinx.mesh import create_unit_square
from ufl import dx, grad, inner
from mpi4py import MPI
from petsc4py import PETSc
from ufl.algorithms.compute_form_data import estimate_total_polynomial_degree
#frequency range definition
f_axis = np.arange(50, 2005, 5)
#Mic position definition
mic = np.array([0.3, 0.3, 0])
p_mic = np.zeros((len(f_axis),1),dtype=complex)
# approximation space polynomial degree
deg = 2
# number of elements in each direction of msh
n_elem = 64
msh = create_unit_square(MPI.COMM_SELF, n_elem, n_elem)
n = ufl.FacetNormal(msh)
# Source amplitude
Q = 1
# Test and trial function space
V = FunctionSpace(msh, ("CG", deg))
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
f = Function(V)
#Source definition position = (Sx,Sy)
Sx = 0.5
Sy = 0.5
alfa = 0.015
delta_tmp = Function(V)
delta_tmp.interpolate(lambda x : 1/(np.abs(alfa)*np.sqrt(np.pi))*np.exp(-(((x[0]-Sx)**2+(x[1]-Sy)**2)/(alfa**2))))
int_delta_tmp = assemble_scalar(form(delta_tmp*dx)) #form() l'ho scoperto per sbaglio. Senza esplode.
delta = delta_tmp/int_delta_tmp
int_delta = assemble_scalar(form(delta*dx))
# fluid definition
c0 = 340
rho_0 = 1.225
omega = Constant(msh, PETSc.ScalarType(1))
k0 = Constant(msh, PETSc.ScalarType(1))
#building the weak form
f = 1j*rho_0*omega*Q*delta
a = inner(grad(u), grad(v)) * dx - k0**2 * inner(u, v) * dx
L = inner(f, v) * dx
#building the problem
uh = Function(V)
uh.name = "u"
problem = LinearProblem(a, L, u=uh, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
# frequency for loop
print(" ")
for nf in range(0,len(f_axis)):
print ("\033[A \033[A")
freq = f_axis[nf]
print("Computing frequency: " + str(freq) + " Hz...")
# Compute solution
omega.value =f_axis[nf]*2*np.pi
k0.value = 2*np.pi*f_axis[nf]/c0
problem.solve()
#Microphone pressure at specified point evaluation
points = np.zeros((3,1))
points[0][0] = mic[0]
points[1][0] = mic[1]
points[2][0] = mic[2]
bb_tree = geometry.BoundingBoxTree(msh, msh.topology.dim)
cells = []
points_on_proc = []
# Find cells whose bounding-box collide with the the points
cell_candidates = geometry.compute_collisions(bb_tree, points.T)
# Choose one of the cells that contains the point
colliding_cells = geometry.compute_colliding_cells(msh, cell_candidates, points.T)
for i, point in enumerate(points.T):
if len(colliding_cells.links(i))>0:
points_on_proc.append(point)
cells.append(colliding_cells.links(i)[0])
points_on_proc = np.array(points_on_proc, dtype=np.float64)
u_values = uh.eval(points_on_proc, cells)
#building of sound pressure vector
p_mic[nf] = u_values
print ("\033[A \033[A")
print("JOB COMPLETED")
#creation of file [frequency Re(p) Im(p)]
Re_p = np.real(p_mic)
Im_p = np.imag(p_mic)
arr = np.column_stack((f_axis,Re_p,Im_p))
print("Saving data to .csv file...")
np.savetxt("p_mic_2nd.csv", arr, delimiter=",")
#spectrum plot
import matplotlib.pyplot as plt
fig = plt.figure()
plt.plot(f_axis, np.log10(np.abs(p_mic)), "k", linewidth=1, label="Sound pressure")
plt.grid(True)
plt.xlabel("f")
plt.legend()
plt.show()
Finally, I want to share the comparison between FEniCSx and Actran, in which the same mesh, elements and polynomial degree have been used (unfortunately I don’t know anything about Actran’s quadrature schemes).
1st order polynomial (TRIA3):
2nd order polynomial (TRIA6):
I hope that this helps someone.