Complex Helmholtz Equation with Neumann Conditions

[stefano.tronci]

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

• Solving complex-valued PDE

Thanks!

[dparsons]

Hi Stefano, what host system are you running on?

[stefano.tronci]

Hi @dparsons , thank you very much for your reply. I am running on Ubuntu 18.04 x86_64. I installed it as in explained at page 7 of this book.

[dparsons]

OK. Ubuntu 18.04 bionic doesn’t support dolfinx, which needs a more recent version of PETSc.

Ubuntu 20.04 focal provides PETSc 3.12, which is known to generate null imaginary components in complex calculations for dolfinx, see complex number tests generate zero solution with PETSc 3.12 (ubuntu focal 20.04 PPA build) · Issue #1612 · FEniCS/dolfinx · GitHub

But since you’re using dolfin not dolfinx, and using explicitly separated real and imaginary components, it would be odd if your zero solution is appearing for the same reason that it happens for dolfinx with petsc 3.12. More likely it’s happening for a different reason.

In a roundabout way what I’m getting at is that you might find it easier to obtain your correct solution using dolfinx, since it supports complex numbers directly (when built with complex number support), no need to split the components. But to use dolfinx successfully you’d need to upgrade to a more recent Ubuntu, 21.04 hirsute or later.

[stefano.tronci]

@dparsons , thank you very much for your insight. Unfortunately at the moment I cannot upgrade my machine to Ubuntu 18.04, but I guess I can give it a spin in some other way (maybe a virtual machine) and see how it goes.

If there are some Ubuntu 18.04 options to try to get this problem working I will be happy to try them.

[conpierce8]

Hi @stefano.tronci,

The reason for the null solution is the use of fenics.ds(i) in your linear form L. You should instead use ds(i) (without “fenics.”). The reason is because you define

@stefano.tronci

ds = fenics.Measure('ds', domain=mesh, subdomain_data=mf)

which is a variable that exists in the namespace of your Python script, while fenics.ds exists within the namespace of the package fenics.

As an aside, your linear and bilinear forms might be missing terms, but I’m not an expert in complex-valued forms so I can’t say for sure. Starting from the strong form in the Elmer link you provided:

P_{,ii} + k^2 P = 0

and introducing the (complex-valued) test function \hat{P}, we can obtain the weak form:

\int_V { \left(P_{,ii} \overline{\hat{P}} + k^2 P \overline{\hat{P}}\right) dx} = 0

where \overline{(\cdot)} is the complex conjugate of (\cdot). Integrating by parts yields:

\int_{\partial V} {P_{,i} n_i \overline{\hat{P}} ds} - \int_V {P_{,i} \overline{\hat{P}_{,i}} dx } + \int_V {k^2 P \overline{\hat{P}} dx } = 0

\int_V {P_{,i} \overline{\hat{P}_{,i}} dx } - \int_V {k^2 P \overline{\hat{P}} dx } = \int_{\partial V} {g \overline{\hat{P}} ds}

Splitting P=P_R + i P_I and \hat{P}=\hat{P}_R + i \hat{P}_I and g = g_R + i g_I gives:

\int_V {P_{,i} \overline{\hat{P}_{,i}} dx } = \int_V { [P_{R,i} \hat{P}_{R,i} + P_{I,i} \hat{P}_{I,i} + i(P_{I,i} \hat{P}_{R,i} - P_{R,i} \hat{P}_{I,i})] dx} \\ \int_V {k^2 P \overline{\hat{P}} dx } = \int_V { k^2 [P_{R} \hat{P}_{R} + P_{I} \hat{P}_{I} + i(P_{I} \hat{P}_{R} - P_{R} \hat{P}_{I})] dx} \\ \int_V {g \overline{\hat{P}} dx } = \int_V { [g_{R} \hat{P}_{R} + g_{I} \hat{P}_{I} + i(g_{I} \hat{P}_{R} - g_{R} \hat{P}_{I})] dx}

So the bilinear form would be:

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 + \
    fenics.inner(fenics.nabla_grad(u_im), fenics.nabla_grad(v_re)) * fenics.dx - k_sq * u_im * v_re * fenics.dx - \
    fenics.inner(fenics.nabla_grad(u_re), fenics.nabla_grad(v_im)) * fenics.dx + k_sq * u_re * v_im * fenics.dx

and your linear form would be:

L = \
    (g_rig_re * v_re + g_rig_im * v_im + g_rig_im * v_re - g_rig_re * v_im) * ds(1) + \
    (g_rig_re * v_re + g_rig_im * v_im + g_rig_im * v_re - g_rig_re * v_im) * ds(2) + \
    (g_rig_re * v_re + g_rig_im * v_im + g_rig_im * v_re - g_rig_re * v_im) * ds(3) + \
    (g_rig_re * v_re + g_rig_im * v_im + g_rig_im * v_re - g_rig_re * v_im) * ds(4) + \
    (g_rig_re * v_re + g_rig_im * v_im + g_rig_im * v_re - g_rig_re * v_im) * ds(5) + \
    (g_in_re * v_re + g_in_im * v_im + g_in_im * v_re - g_in_re * v_im) * ds(6)

Hopefully someone from the community will correct me if this is incorrect.

Regards,
Connor

[stefano.tronci]

Hi @conpierce8 , thank you very much for your help! I see I did a rookie mistake with the ds, classic. By correcting the ds definition the imaginary part of the solution matches that of Elmer (qualitatively, the values are off). The real part is still zero.

If, on top of that, I replace my bilinear and linear forms with yours I still get a wrong solution: identical real and imaginary parts. Maybe that’s the right direction though. A bilinear form for this problem was provided here. I am not sure whether that is correct, but I guess I will try to make sure I am implementing that and see how it goes.

[conpierce8]

Hi @stefano.tronci,

The zero solution for the real part may be correct, since you have a lossless material and the Neumann condition is purely imaginary. Does your Elmer solution also assume a lossless material? I note that the tutorial mentions that it is possible to include a damping factor D.

As for the weak form, it seems that the accepted solution in the thread you linked gives the same weak form that I derived. The post by @BillS is stated to be for a lossless medium; perhaps these two weak forms become equivalent when there is no damping?

[stefano.tronci]

Hi @conpierce8 ,

Good point. I can confirm D is set to 0 in my Elmer simulation, so I reckon the two solvers are dealing with the same equation.

[conpierce8]

Upon looking closer at your Elmer results, it appears that the imaginary part is negligible O(10^{-19}) compared to the real part O(10^1). Did you apply a real-valued Neumann BC in Elmer?

[BillS]

Yes, for the lossless problem, the two weak forms are identical. This is because the real and imaginary solutions to the PDE are identical for the lossless medium and the cross terms will cancel out. Real and imaginary solutions are decoupled. Occam’s Razor is used to cut away the unneeded terms that only enlarge the matrix equation. You only need to have the cross terms when losses are present in the volume or surface element under consideration.

[stefano.tronci]

@conpierce8

Ah dang, yes I did apply a real-valued Neumann BC for that plot. I was playing with it and I forgot to reset it to the same value as in the Fenics code. Below is the Elmer solution with matching imaginary Neumann BC:

[stefano.tronci]

Not only that, but I did not apply the Neumann condition to the same boundary in the two simulations… This is the Fenics solution (complete weak form with cross terms) when the influx boundary condition is applied to the same boundary as in the Elmer simulation:

Seems like things are looking pretty good. Results are very similar. Thank you for bearing with me guys, I think this is solved!