Hi!
I am working with a slightly more general version of the Helmholtz equation in Acoustics. My equation just includes two material damping parameters, \eta and \alpha. I derived the weak form here. I have implemented a FEniCS solver for it but it provides the wrong result. The details of my problem are below.
To verify that things work I have created a simple domain. It is a sphere of 0.1 meters radius. I cut off a sphere of 5 mm radius from the center of it. A slice of the mesh (Salome) is shown below.
We can call the inner surface (5 mm radius) surface 1. We can call the outer surface (0.1 m radius) surface 0. I want to apply an uniform velocity to surface 1. This is a Neumann boundary condition which will make the inner surface a pulsating sphere. Then I want to apply a matched impedance to surface 0 so that there will be no reflections. This is a Robin boundary condition.
My weak form is as follows:
- \gamma \int_{\Omega} \nabla u \cdot \nabla \phi^\star d\mathbf{x} + \beta \int_{\Omega} u \phi^{\star} d\mathbf{x} + \frac{j \omega \rho \gamma}{z \left( 0.1 \text{ m} \right)} \int_{\partial \Omega_{0}} u \phi^\star ds = - j \omega \rho \gamma U \int_{\partial \Omega_{1}} \phi^\star ds
Where:
- \Omega is the domain, \partial\Omega_{0} and \partial\Omega_{1} are the outer and inner surfaces respectively;
- u: \Omega \rightarrow \mathbb{C} is the trial function;
- \phi: \Omega \rightarrow \mathbb{C} is the test function;
- \gamma \doteq 1 + j k \frac{\eta}{c}, with \eta a damping parameter (set to 0 for this simulation);
- \beta \doteq k^2 - jk \frac{\alpha}{c}, with \alpha a damping parameter (set to 0 for this simulation);
- k = \frac{\omega}{c} = \frac{2 \pi \nu}{c}, with \nu the frequency of the simulation (1000 Hz) and c the speed of sound in air (343 m/s);
- z is the uniform specific impedance of outgoing spherical waves: z\left(r\right)=\rho c \frac{jkr}{1+jkr}, with \rho the density of air (1.205 kg/m^3);
- U is the uniform normal velocity of the inner surface, set to \frac{\sqrt{2}}{2} + j\frac{\sqrt{2}}{2} m/s for this simulation.
- j is the imaginary unit. \star denotes complex conjugation.
z and U pop out of the integrals due to being uniform. Since \eta = \alpha = 0 I expect the solution to be that of the lossless pulsating sphere:
u\left(r\right) = z\left(a \right)\frac{aU}{r}\exp\left(-jk \left( r-a\right) \right)
With a the radius of the source (5 mm). However, my implementation below is quite wrong (sorry for it being 222 lines of code, this is the most I could squeeze it):
import pathlib
import meshio
import tempfile
import numpy as np
import fenics
import matplotlib.pyplot as plt
# This is a 3D domain which is a sphere of 0.1 m radius.
# A smaller sphere of 5 mm radius is cutout from its center.
# We will apply a uniform velocity (Neumann) boundary condition to the inner surface.
# We will apply an impedance (Robin) boundary condition to the outer surface.
# Since the mesh is in mm, we convert to m with `coordinate_scaling`.
mesh_path = pathlib.Path('spherecut.med')
coordinate_scaling = 1.0e-3
# We read te MED file with meshio
raw_mesh_data = meshio.read(mesh_path)
# MED cell tags are negative, but FEniCS doesn't like it. We convert to positive by applying this offset.
cell_tag_offset = np.abs(np.min(list(raw_mesh_data.cell_tags.keys())))
# This list will contain FEniCS meshes, in this order: line (1D), triangle (2D) and tetra (3D).
# To make the FEniCS meshes we crate meshio meshes for each entity, save it to a temporary location as xdmf, and then
# we read it back.
mesh = fenics.Mesh()
mvc_3d = fenics.MeshValueCollection('size_t', mesh, 3)
mvc_2d = fenics.MeshValueCollection('size_t', mesh, 2)
mvc_1d = fenics.MeshValueCollection('size_t', mesh, 1)
with tempfile.TemporaryDirectory() as xdmf_dir:
target_xdmf_dir = pathlib.Path(xdmf_dir)
for entity in ['tetra', 'triangle', 'line']:
cells = raw_mesh_data.get_cells_type(entity)
cell_data = raw_mesh_data.get_cell_data('cell_tags', entity)
points = raw_mesh_data.points # 3D Mesh: No pruning needed.
edited_mesh_data = meshio.Mesh(
points=points * coordinate_scaling,
cells={entity: cells},
cell_data={'marker': [cell_data + cell_tag_offset]}
)
tmp_path = target_xdmf_dir.joinpath(entity + '_mesh.xdmf')
meshio.write(tmp_path.absolute().as_posix(), edited_mesh_data)
with fenics.XDMFFile(tmp_path.absolute().as_posix()) as xdmf_file:
if entity == 'tetra':
xdmf_file.read(mesh)
xdmf_file.read(mvc_3d, 'marker')
if entity == 'triangle':
xdmf_file.read(mvc_2d, 'marker')
if entity == 'line':
xdmf_file.read(mvc_1d, 'marker')
markers = fenics.MeshFunction('size_t', mesh, mvc_3d)
sub_domains_2d = fenics.MeshFunction('size_t', mesh, mvc_2d)
sub_domains_1d = fenics.MeshFunction('size_t', mesh, mvc_1d)
# We have the mesh, now for the simulation. The simulation parameters are below.
frequency = 1000.0
speed_of_sound = 343.0
density = 1.205
eta = 0.0
alpha = 0.0
source_velocity = (1.0 + 1.0j) / np.sqrt(2.0)
# We now define the function space. The PDE is complex, so we use mixed space.
base_element = fenics.FiniteElement(
family='Lagrange',
cell=mesh.ufl_cell(),
degree=1
)
# element = fenics.MixedElement([base_element, base_element])
# f_space = fenics.FunctionSpace(mesh, element)
f_space = fenics.FunctionSpace(mesh, base_element * base_element)
# Now we can get the trial and test functions (both parts)
trial_re, trial_im = fenics.TrialFunction(V=f_space)
test_re, test_im = fenics.TestFunction(V=f_space)
# We compute a bunch of useful quantities here.
angular_frequency = 2.0 * np.pi * frequency
wave_number = angular_frequency / speed_of_sound
gamma = 1.0 + 1.0j * wave_number * eta / speed_of_sound
beta = wave_number ** 2 - 1.0j * alpha / speed_of_sound
g_flux = 1.0j * angular_frequency * density * gamma
# Nice. Let's start with the first two terms of the bilinear form.
# This is - \gamma \int_{\Omega} \nabla u \cdot \nabla \phi^\star d\mathbf{x}
# (all expanded in real and imaginary parts)
u_gd_re = fenics.grad(trial_re)
u_gd_im = fenics.grad(trial_im)
p_gd_re = fenics.grad(test_re)
p_gd_im = fenics.grad(test_im)
md_dot_re = fenics.dot(u_gd_re, p_gd_re) + fenics.dot(u_gd_im, p_gd_im)
md_dot_im = fenics.dot(u_gd_im, p_gd_re) - fenics.dot(u_gd_re, p_gd_im)
g_md_dot_re = fenics.Constant(-gamma.real)
g_md_dot_im = fenics.Constant(-gamma.imag)
md_dot_term_re = g_md_dot_re * md_dot_re - g_md_dot_im * md_dot_im
md_dot_term_im = g_md_dot_re * md_dot_im + g_md_dot_im * md_dot_re
# This is \beta \int_{\Omega} u \phi^{\star} d\mathbf{x}
# (all expanded in real and imaginary parts)
md_prod_re = trial_re * test_re + trial_im * test_im
md_prod_im = trial_im * test_re - trial_re * test_im
g_md_prod_re = fenics.Constant(beta.real)
g_md_prod_im = fenics.Constant(beta.imag)
md_prod_term_re = g_md_prod_re * md_prod_re - g_md_prod_im * md_prod_im
md_prod_term_im = g_md_prod_re * md_prod_im + g_md_prod_im * md_prod_re
# Now that we have all the bits we sum them together.
bilinear_form = (
md_dot_term_re * fenics.dx +
md_dot_term_im * fenics.dx +
md_prod_term_re * fenics.dx +
md_prod_term_im * fenics.dx
)
# For the boundary conditions, we produce the measure. According to `raw_mesh_data.cell_tags` the outer layer is marked
# with `-7`, the inner with `-6`. We need to add the offset (`7`). So the outer is marked with `0` and the inner with
# `1`.
ds_c = fenics.Measure('ds', domain=mesh, subdomain_data=sub_domains_2d)
# `ds_c(0)` should return a measure over the outer shell.
# `ds_c(1)` should return a measure over the inner shell.
# We also compute the radii of the inner and outer shells straight from the mesh data.
outer_radius = np.max(np.sqrt(np.sum((coordinate_scaling * raw_mesh_data.points) ** 2, axis=1)))
inner_radius = np.min(np.sqrt(np.sum((coordinate_scaling * raw_mesh_data.points) ** 2, axis=1)))
# For the outer shell, we apply an impedance that is matched to that of a spherical wave (no reflections).
# This is a Robin boundary condition.
# This adds this term to the bilinear form:
# j \omega \rho \gamma \int_{\partial \Omega_{0}} \frac{u}{z} \phi^\star ds
jkr = 1.0j * wave_number * outer_radius
impedance = density * speed_of_sound * jkr / (1.0 + jkr)
g_rb = g_flux / impedance
rb_prod_re = trial_re * test_re + trial_im * test_im
rb_prod_im = trial_im * test_re - trial_re * test_im
g_rb_re = fenics.Constant(g_rb.real)
g_rb_im = fenics.Constant(g_rb.imag)
rb_term_re = g_rb_re * rb_prod_re - g_rb_im * rb_prod_im
rb_term_im = g_rb_re * rb_prod_im + g_rb_im * rb_prod_re
bilinear_form += rb_term_re * ds_c(0) + rb_term_im * ds_c(0)
# For the inner shell we apply a normal velocity.
# This is a Neumann boundary condition.
# This adds this term to the linear form:
# - j \omega \rho \gamma \int_{\partial \Omega_{1}} w \phi^\star ds
g_nm = -1.0 * g_flux * source_velocity
g_nm_re = fenics.Constant(g_nm.real)
g_nm_im = fenics.Constant(g_nm.imag)
nm_prod_re = g_nm_re * test_re + g_nm_im * test_im
nm_prod_im = g_nm_im * test_re - g_nm_re * test_im
linear_form = nm_prod_re * ds_c(1) + nm_prod_im * ds_c(1)
# Cool, we are ready to solve!
u = fenics.Function(f_space)
fenics.solve(
bilinear_form == linear_form,
u,
solver_parameters={
'linear_solver': 'umfpack',
# 'linear_solver': 'bicgstab',
# 'preconditioner': 'ilu',
},
)
# To check the results over the mesh nodes I like to write to VTU and read back.
# (it is easier than using the `u_re` and `u_im` objects).
u_re, u_im = u.split(deepcopy=True)
u_re.rename('re', 're')
u_im.rename('im', 'im')
vtk_file_1 = fenics.File('test_re.pvd')
vtk_file_1 << u_re
vtk_file_2 = fenics.File('test_im.pvd')
vtk_file_2 << u_im
u_re_vtu = meshio.read('test_re000000.vtu')
u_im_vtu = meshio.read('test_im000000.vtu')
nodes = u_re_vtu.points
u_re_at_nodes = u_re_vtu.point_data['re']
u_im_at_nodes = u_im_vtu.point_data['im']
u = u_re_at_nodes + 1j * u_im_at_nodes
# We compare with the exact lossless solution.
radii = np.sqrt(np.sum(nodes**2, axis=1))
jka = 1.0j * wave_number * inner_radius
z_a = density * speed_of_sound * jka / (1.0 + jka)
u_exact = z_a * inner_radius * source_velocity * np.exp(-1j * wave_number * (radii - inner_radius)) / radii
idx = np.argsort(radii)
error = u / u_exact
fig = plt.figure('Error Plot')
axs = fig.subplots(2, 1, sharex=True)
axs[0].grid(True, which='both')
axs[1].grid(True, which='both')
axs[1].set_xlabel('r [m]')
axs[0].set_ylabel('Magnitude [dB]')
axs[1].set_ylabel('Phase [rad]')
axs[0].plot(radii[idx], 20 * np.log10(np.abs(error[idx])))
axs[1].plot(radii[idx], np.angle(error[idx]))
plt.show()
fig = plt.figure('Solution Plot')
axs = fig.subplots(2, 1, sharex=True)
axs[0].grid(True, which='both')
axs[1].grid(True, which='both')
axs[1].set_xlabel('r [m]')
axs[0].set_ylabel('SPL [dBSPL]')
axs[1].set_ylabel('Phase [rad]')
axs[0].plot(radii[idx], 20 * np.log10(np.abs(u_exact[idx]) / (20e-6 * np.sqrt(2))))
axs[0].plot(radii[idx], 20 * np.log10(np.abs(u[idx]) / (20e-6 * np.sqrt(2))))
fig.legend(['Exact', 'FEniCS'])
axs[1].plot(radii[idx], np.angle(u_exact[idx]))
axs[1].plot(radii[idx], np.angle(u[idx]))
plt.show()
The mesh to run this code can be downloaded here.
Error plot in magnitude and phase:
Comparison with exact solution:
I have been trying to figure out what I did wrong for a couple of days without success, but I still feel like I am maybe missing something obvious.
It feels like there are reflections coming from the outer shell, so probably something dodgy with my impedance boundary condition. Any idea?
