In time: $$ \frac{1}{c^2} \frac{\partial^2 p}{\partial t^2} - \left ( \frac{\partial^2 p}{\partial x^2} + \frac{\partial^2 p }{\partial y^2} \right ) = s(x, y, t) $$
Assuming $p(x, y, t) = \Re ( P(x, y) e^{j \omega t})$ (this notation is coherent with the Fourier transform notation used habitually, where $s(t) = \int_{-\infty}^{\infty} \hat{s}(\omega) e ^{j \omega t} d\omega$):
$$ \frac{\omega^2}{c^2} P(x, y) + \frac{\partial^2 P}{\partial x^2} + \frac{\partial^2 P}{\partial y^2} = S(x, y) $$
Let's try and setup the Helmholtz equation in 2D.
BC : zero field
import numpy as np
import holoviews as hv
hv.extension('bokeh')
nx, ny = 5, 5
x = np.linspace(0, 1, num=nx)
y = np.linspace(0, 1, num=ny)
dx = (x[-1] - x[0]) / nx
dy = (y[-1] - y[0]) / ny
Let's plot the coordinates of each point.
points = []
for i in range(nx):
for j in range(ny):
xx = x[i]
yy = y[j]
points.append((xx, yy))
labels = hv.Labels({('x', 'y'): points, 'Label': [chr(65+i) for i in range(len(points))]}).opts(xlim=(-.1, 1.1), ylim=(-.1, 1.1))
labels
Let's now write the main loop for assembling the matrix that will hold the system.
def ij_to_row_coords(i, j, nx, ny):
return ny * i + j
def row_to_ij_coords(row, nx, ny):
i = row // ny
j = row - i * ny
return i, j
omega = 2 * np.pi * 800
c = 1500
k = omega / c
N = nx * ny
F = np.zeros((N, N), dtype=complex)
for i in range(nx):
for j in range(ny):
if i in [0, nx - 1] or j in [0, ny - 1]:
if i == 0:
# left
ij = ij_to_row_coords(i, j, nx, ny)
ip1j = ij_to_row_coords(i+1, j, nx, ny)
F[ij, ij] += 1/dx - 1j*k
F[ij, ip1j] += -1/dx
if i == nx - 1:
# right
ij = ij_to_row_coords(i, j, nx, ny)
im1j = ij_to_row_coords(i-1, j, nx, ny)
F[ij, ij] += 1/dx - 1j*k
F[ij, im1j] += -1/dx
if j == 0:
# bottom
ij = ij_to_row_coords(i, j, nx, ny)
ijp1 = ij_to_row_coords(i, j+1, nx, ny)
F[ij, ij] += 1/dx - 1j * k
F[ij, ijp1] += -1/dx
if j == ny - 1:
# top
ij = ij_to_row_coords(i, j, nx, ny)
ijm1 = ij_to_row_coords(i, j-1, nx, ny)
F[ij, ij] += 1/dx - 1j * k
F[ij, ijm1] += -1/dx
else:
# interior
ij = ij_to_row_coords(i, j, nx, ny)
ip1j = ij_to_row_coords(i+1, j, nx, ny)
im1j = ij_to_row_coords(i-1, j, nx, ny)
ijp1 = ij_to_row_coords(i, j+1, nx, ny)
ijm1 = ij_to_row_coords(i, j-1, nx, ny)
F[ij, ij] += k**2 - 2/dx**2 - 2/dy**2
F[ij, ip1j] += 1/dx**2
F[ij, im1j] += 1/dx**2
F[ij, ijp1] += 1/dy**2
F[ij, ijm1] += 1/dy**2
Let's visualize F:
F.min(), F.max()
hv.Image(np.real(F)).opts(colorbar=True, width=350) + hv.Image(np.imag(F)).opts(colorbar=True, width=350)
We can also visualize non-zero coordinates as a function of their location.
def show_nonzero(line):
nonzeros = F[line].nonzero()
center_ij = row_to_ij_coords(line, nx, ny)
center_point = hv.Points(np.c_[(x[center_ij[0]], y[center_ij[1]])]).opts(size=20)
center_label = hv.Text(x[center_ij[0]], y[center_ij[1]], str(center_ij))
nonzero_ijs = [row_to_ij_coords(row, nx, ny) for row in nonzeros[0]]
nonzero_points = hv.Points(np.c_[[(x[i], y[j]) for i,j in nonzero_ijs]]).opts(size=10)
return center_point * nonzero_points * center_label
labels * hv.DynamicMap(show_nonzero, kdims='line').redim.range(line=(0, nx * ny - 1))
The equation to solve is $F \cdot u = s$.
s = np.zeros(F.shape[0])
s[ij_to_row_coords(1, 1, nx, ny)] = -1
s
Let's now try to solve this system.
np.linalg.solve(F, s)
Let's move on and simulate the realistic case shown in the paper. The grid is now a 400 × 300 grid with ∆x = ∆z = 0.11m. The paper explains that the discretization has been chosen to verify ∆x < λ/12.
nx = 400
ny = 300
x = np.linspace(0, 44, num=nx)
y = np.linspace(0, 33, num=ny)
dx = (x[-1] - x[0]) / nx
dy = (y[-1] - y[0]) / ny
omega = 2 * np.pi * 800
c = 1500
k = omega / c
nx, ny, dx, dy, k, x.min(), x.max(), y.min(), y.max()
This time, we can write a function to assemble the matrix.
from scipy.sparse.linalg import spsolve
from scipy.sparse import lil_matrix, csc_matrix
from tqdm import tqdm_notebook
def assemble_F_matrix(nx, ny, k, dx):
"""Assembles the F matrix to solve the Helmholtz equation.
Assumes dx = dy.
"""
N = nx * ny
F = lil_matrix((N, N), dtype=complex)
for i in tqdm_notebook(range(nx)):
for j in range(ny):
if i in [0, nx - 1] or j in [0, ny - 1]:
if i == 0:
# left
ij = ij_to_row_coords(i, j, nx, ny)
ip1j = ij_to_row_coords(i+1, j, nx, ny)
F[ij, ij] += 1/dx - 1j*k
F[ij, ip1j] += -1/dx
if i == nx - 1:
# right
ij = ij_to_row_coords(i, j, nx, ny)
im1j = ij_to_row_coords(i-1, j, nx, ny)
F[ij, ij] += 1/dx - 1j*k
F[ij, im1j] += -1/dx
if j == 0:
# bottom
ij = ij_to_row_coords(i, j, nx, ny)
ijp1 = ij_to_row_coords(i, j+1, nx, ny)
F[ij, ij] += 1/dx - 1j * k
F[ij, ijp1] += -1/dx
if j == ny - 1:
# top
ij = ij_to_row_coords(i, j, nx, ny)
ijm1 = ij_to_row_coords(i, j-1, nx, ny)
F[ij, ij] += 1/dx - 1j * k
F[ij, ijm1] += -1/dx
else:
# interior
ij = ij_to_row_coords(i, j, nx, ny)
ip1j = ij_to_row_coords(i+1, j, nx, ny)
im1j = ij_to_row_coords(i-1, j, nx, ny)
ijp1 = ij_to_row_coords(i, j+1, nx, ny)
ijm1 = ij_to_row_coords(i, j-1, nx, ny)
F[ij, ij] += k**2 - 2/dx**2 - 2/dx**2
F[ij, ip1j] += 1/dx**2
F[ij, im1j] += 1/dx**2
F[ij, ijp1] += 1/dx**2
F[ij, ijm1] += 1/dx**2
return F
Let's now see if we can model a simple point source (which is a line source in 3D).
F = assemble_F_matrix(nx, ny, k, dx)
i_source = np.argmin(np.abs(x - 11.))
j_source = np.argmin(np.abs(y - 11.))
s = np.zeros(F.shape[0], dtype=complex)
s[ij_to_row_coords(i_source, j_source, nx, ny)] = -1.
solution = spsolve(csc_matrix(F), s)
We can now unpack the solution to display it.
abs_sol = np.abs(solution).reshape((nx, ny)).T
real_sol = np.real(solution).reshape((nx, ny)).T
imag_sol = np.imag(solution).reshape((nx, ny)).T
zmin, zmax = np.abs(solution).min(), np.abs(solution).max()
zval = zmax / 1
img_opts = hv.opts.Image(colorbar=False, width=400, height=300)
(hv.Image(abs_sol, bounds=(x.min(), y.min(), x.max(), y.max()), label='abs').opts(img_opts).redim.range(z=(-zval, zval)) + \
hv.Image(real_sol, bounds=(x.min(), y.min(), x.max(), y.max()), label='real').opts(img_opts).redim.range(z=(-zval, zval)) + \
hv.Image(imag_sol, bounds=(x.min(), y.min(), x.max(), y.max()), label='imag').opts(img_opts).redim.range(z=(-zval, zval))).cols(2)
What does this become if we change the frequency?
omega = 2 * np.pi * 400
k = omega / c
F = assemble_F_matrix(nx, ny, k, dx)
solution = spsolve(csc_matrix(F), s)
abs_sol = np.abs(solution).reshape((nx, ny)).T
real_sol = np.real(solution).reshape((nx, ny)).T
imag_sol = np.imag(solution).reshape((nx, ny)).T
zmin, zmax = np.abs(solution).min(), np.abs(solution).max()
zval = zmax / 10
img_opts = hv.opts.Image(colorbar=False, width=400, height=300)
(hv.Image(abs_sol, bounds=(x.min(), y.min(), x.max(), y.max()), label='abs').opts(img_opts).redim.range(z=(-zval, zval)) + \
hv.Image(real_sol, bounds=(x.min(), y.min(), x.max(), y.max()), label='real').opts(img_opts).redim.range(z=(-zval, zval)) + \
hv.Image(imag_sol, bounds=(x.min(), y.min(), x.max(), y.max()), label='imag').opts(img_opts).redim.range(z=(-zval, zval))).cols(2)
As expected, we get a lower frequency solution.
We can also experiment with a different source. Let's try a piston like source.
i_source = np.argmin(np.abs(x - 1.))
s = np.zeros(F.shape[0], dtype=complex)
for j_source in range(1, ny - 1):
s[ij_to_row_coords(i_source, j_source, nx, ny)] = -1.
solution = spsolve(csc_matrix(F), s)
abs_sol = np.abs(solution).reshape((nx, ny)).T
real_sol = np.real(solution).reshape((nx, ny)).T
imag_sol = np.imag(solution).reshape((nx, ny)).T
zmin, zmax = np.abs(solution).min(), np.abs(solution).max()
zval = zmax
img_opts = hv.opts.Image(colorbar=False, width=400, height=300)
(hv.Image(abs_sol, bounds=(x.min(), y.min(), x.max(), y.max()), label='abs').opts(img_opts).redim.range(z=(-zval, zval)) + \
hv.Image(real_sol, bounds=(x.min(), y.min(), x.max(), y.max()), label='real').opts(img_opts).redim.range(z=(-zval, zval)) + \
hv.Image(imag_sol, bounds=(x.min(), y.min(), x.max(), y.max()), label='imag').opts(img_opts).redim.range(z=(-zval, zval))).cols(2)
by changing the size of the source with respect to the wavelength, we can see a difference in source size.
Above, the lenght of the source was 33 meters at 400 Hz which yields, in terms of wavelengths:
D_source = ny * dx
wavelength = 2 * np.pi * c / omega
D_source / wavelength
What if we make the source equal to 0.5 times the wavelength?
0.5 * wavelength / dx
i_source = np.argmin(np.abs(x - 1.))
s = np.zeros(F.shape[0], dtype=complex)
for j_source in range(int(ny/2 - 8), int(ny/2 + 8)):
s[ij_to_row_coords(i_source, j_source, nx, ny)] = -1.
solution = spsolve(csc_matrix(F), s)
abs_sol = np.abs(solution).reshape((nx, ny)).T
real_sol = np.real(solution).reshape((nx, ny)).T
imag_sol = np.imag(solution).reshape((nx, ny)).T
zmin, zmax = np.abs(solution).min(), np.abs(solution).max()
zval = zmax
img_opts = hv.opts.Image(colorbar=False, width=400, height=300)
(hv.Image(abs_sol, bounds=(x.min(), y.min(), x.max(), y.max()), label='abs').opts(img_opts).redim.range(z=(-zval, zval)) + \
hv.Image(real_sol, bounds=(x.min(), y.min(), x.max(), y.max()), label='real').opts(img_opts).redim.range(z=(-zval, zval)) + \
hv.Image(imag_sol, bounds=(x.min(), y.min(), x.max(), y.max()), label='imag').opts(img_opts).redim.range(z=(-zval, zval))).cols(2)
Interestingly we get a radiation pattern that ressembles that of a point source.
Finally, we can get an animation of the complex radiated field by an intermediate source of size 4 wavelengths.
4 * wavelength / dx
i_source = np.argmin(np.abs(x - 1.))
s = np.zeros(F.shape[0], dtype=complex)
for j_source in range(int(ny/2 - 70), int(ny/2 + 70)):
s[ij_to_row_coords(i_source, j_source, nx, ny)] = -1.
solution = spsolve(csc_matrix(F), s)
abs_sol = np.abs(solution).reshape((nx, ny)).T
real_sol = np.real(solution).reshape((nx, ny)).T
imag_sol = np.imag(solution).reshape((nx, ny)).T
zmin, zmax = np.abs(solution).min(), np.abs(solution).max()
zval = zmax
img_opts = hv.opts.Image(colorbar=False, width=400, height=300)
(hv.Image(abs_sol, bounds=(x.min(), y.min(), x.max(), y.max()), label='abs').opts(img_opts).redim.range(z=(-zval, zval)) + \
hv.Image(real_sol, bounds=(x.min(), y.min(), x.max(), y.max()), label='real').opts(img_opts).redim.range(z=(-zval, zval)) + \
hv.Image(imag_sol, bounds=(x.min(), y.min(), x.max(), y.max()), label='imag').opts(img_opts).redim.range(z=(-zval, zval))).cols(2)
Let's synthesize the associated time varying field.
snapshots = {}
timesteps = 30
for timestep in range(timesteps):
t = timestep / timesteps * 2 * np.pi / omega
field = (real_sol + 1j * imag_sol) * np.exp(-1j * omega * t)
snapshots[timestep] = hv.Image(np.real(field), bounds=(x.min(), y.min(), x.max(), y.max()), label='abs').opts(img_opts).redim.range(z=(-zval, zval))
%%output holomap='scrubber'
hv.HoloMap(snapshots)