Hello!

As the title says, I have been trying to solve the Helmholtz equation in a box with a different Neumann condition on each wall. The box represents a room in which harmonic acoustic wave propagation is happening. I used a number of other posts as a guide, but no success. See my code below.

```
import numpy as np
import fenics
# Room Acoustics: Helmholtz Equation in a room
f = 57.17 # Frequency of the Simulation, Hz
c = 343 # Speed of sound in air, m/s
rho = 1.205 # Density of air, kg/m^3
v_n = 10 # Velocity normal to boundary, for inflow Neumann condition, m/s
lx = 4 # Room size x, m
ly = 5 # Room size y, m
lz = 3 # Room size z, m
tol = 1e-10 # Tolerance for boundary condition definitions
# Computing some useful stuff
w = 2 * np.pi * f # Angular frequency, rad/s
k = w / c # Wave number, rad/m
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))
Nz = np.intp(np.ceil(lz / s))
# Making a mesh, then an element. Then, Mixed Space so that we can solve for real and imaginary parts.
mesh = fenics.BoxMesh(
fenics.Point(0, 0, 0),
fenics.Point(lx, ly, lz),
Nx,
Ny,
Nz
)
P = fenics.FiniteElement("Lagrange", mesh.ufl_cell(), 2)
V = fenics.FunctionSpace(mesh, P * P)
# I also tried this, but same results:
# element = fenics.MixedElement([P, P])
# V = fenics.FunctionSpace(mesh, element)
# Define variational problem
u_re, u_im = fenics.TrialFunctions(V)
v_re, v_im = fenics.TestFunctions(V)
k_sq = fenics.Constant(k**2)
# As far as I understand this is the correct bilinear form. Right?
a = \
fenics.inner(fenics.nabla_grad(u_re), fenics.nabla_grad(v_re)) * fenics.dx - k_sq * v_re * u_re * fenics.dx + \
fenics.inner(fenics.nabla_grad(u_im), fenics.nabla_grad(v_im)) * fenics.dx - k_sq * v_im * u_im * fenics.dx
# Now I want to be able to use a different Neumann condition on each wall (even though g will be 0 on all walls aside
# for one...)
mf = fenics.MeshFunction("size_t", mesh, 2)
mf.set_all(0)
class BX0(fenics.SubDomain):
def inside(self, x, on_boundary):
return on_boundary and fenics.near(x[0], 0, tol)
class BXL(fenics.SubDomain):
def inside(self, x, on_boundary):
return on_boundary and fenics.near(x[0], lx, tol)
class BY0(fenics.SubDomain):
def inside(self, x, on_boundary):
return on_boundary and fenics.near(x[1], 0, tol)
class BYL(fenics.SubDomain):
def inside(self, x, on_boundary):
return on_boundary and fenics.near(x[1], ly, tol)
class BZ0(fenics.SubDomain):
def inside(self, x, on_boundary):
return on_boundary and fenics.near(x[2], 0, tol)
class BZL(fenics.SubDomain):
def inside(self, x, on_boundary):
return on_boundary and fenics.near(x[2], lz, tol)
bx0 = BX0()
bx0.mark(mf, 1)
bxl = BXL()
bxl.mark(mf, 2)
by0 = BY0()
by0.mark(mf, 3)
byl = BYL()
byl.mark(mf, 4)
bz0 = BZ0()
bz0.mark(mf, 5)
bzl = BZL()
bzl.mark(mf, 6)
# The flux is 0 for all rigid walls
flux_rig = 0
g_rig_re = fenics.Constant(np.real(flux_rig))
g_rig_im = fenics.Constant(np.imag(flux_rig))
# The flux is this for the active walls
flux_in = 1j * w * rho * v_n
g_in_re = fenics.Constant(np.real(flux_in))
g_in_im = fenics.Constant(np.imag(flux_in))
ds = fenics.Measure('ds', domain=mesh, subdomain_data=mf)
# This produces a solution which has null imaginary part and null real part. Why?
L = \
g_rig_re * v_re * fenics.ds(1) + g_rig_im * v_im * fenics.ds(1) + \
g_rig_re * v_re * fenics.ds(2) + g_rig_im * v_im * fenics.ds(2) + \
g_rig_re * v_re * fenics.ds(3) + g_rig_im * v_im * fenics.ds(3) + \
g_rig_re * v_re * fenics.ds(4) + g_rig_im * v_im * fenics.ds(4) + \
g_rig_re * v_re * fenics.ds(5) + g_rig_im * v_im * fenics.ds(5) + \
g_in_re * v_re * fenics.ds(6) + g_in_im * v_im * fenics.ds(6)
# Tried this. Of course it makes all walls having the influx boundary condition. The real part makes kind of sense,
# but the imaginary part is null. Why?
# L = g_in_re * v_re * fenics.ds + g_in_im * v_im * fenics.ds
# Compute solution
u = fenics.Function(V)
fenics.solve(a == L, u)
# Save solution to file in VTK format
u_1, u_2 = u.split()
u_1.rename('re', 're')
u_2.rename('im', 'im')
vtk_file_1 = fenics.File('helmholtz/solution_re.pvd')
vtk_file_1 << u_1
vtk_file_2 = fenics.File('helmholtz/solution_im.pvd')
vtk_file_2 << u_2
```

For reference, I am pretty much attempting to implement the Helmoltz model as in Elmer. This simulation, when setup in Elmer, produces the following results:

However my Fenics code will make a null solution. I think I am missing something vary basic, perhaps in how `ds`

is defined. But I am not completely sure about the variational formulation either. Any idea?

Here the main post I used as a guide (there are other, but I can only post 2 links…):

Thanks!