2D Acoustic simulation: monopole source in a square room - FEniCS vs. Actran

1

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:

MESH_PALLOZZ
MESH_PALLOZZ712×645 130 KB

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):

tria3
tria3990×503 165 KB

2nd order polynomial (TRIA6):

tria6
tria6998×506 166 KB

I hope that this helps someone.

8 Likes
2

If anyone wants a deeper explanation of the steps of the code, you can find it in my blog about computational acoustics:

https://undaproject.com/

4 Likes
3

Thank you very much for your work. Always a pleasure to find the exact case we are looking for, and you blog is full of them.
The difference between the curves in the first plot comes from a bad meshing procedure. In acoustic, most of the time using quadratic element is needed (It’s not always the case, but in acoustic a coarser mesh with high order of element will lead to more accurate result).

Thanks again for your work !

4

First, thank you for your kind words.

You are almost right! In acoustics, it is well known that linear elements produce waves that are “too stiff”.

Here I intentionally used the same (poor) mesh, so in theory, even if the result is inaccurate, under identical conditions all solvers should produce nearly the same output.

However, Actran by default employs an implementation that reduces phase (dispersion) error for first-order elements. As a result, on the same (poor) mesh Actran appears much less “stiff”, though still inaccurate compared to second-order elements.

In conclusion: yes, you are absolutely right, always use quadratic elements!

1 Like