9b737b8GPU implementation of fmm3dbie's interior Helmholtz impedance solver, with a showcase room problemJeremy Magland 1"""Plot the room field saved by run.py (python plot.py results/<case>-<backend>): slices of Re(u) and level in dB,
2plus an animation of the real pressure field Re(u exp(-i omega t)) over one period."""
4import argparse
5import os
7import matplotlib
9matplotlib.use("Agg")
10import matplotlib.pyplot as plt
11import numpy as np
12from matplotlib.animation import FuncAnimation, PillowWriter
14SLICES = [("xy", "x", "y", 2), ("xz", "x", "z", 1), ("yz", "y", "z", 0)]
17def masked(a):
18 return np.ma.masked_invalid(a)
21def main():
22 ap = argparse.ArgumentParser()
23 ap.add_argument("out", help="a results directory written by run.py")
24 ap.add_argument("--frames", type=int, default=32)
25 args = ap.parse_args()
27 d = np.load(os.path.join(args.out, "field.npz"))
28 g, x0, k = d["grid"], d["x0"], float(d["k"])
29 ext = [g[0], g[-1], g[0], g[-1]]
30 us = {n: d[f"u_{n}"] for n, *_ in SLICES}
32 # Color limits from the field more than half a wavelength from the source
33 A, B = np.meshgrid(g, g, indexing="xy")
34 far = []
35 for n, lx, ly, fixed in SLICES:
36 r = np.hypot(A - x0["xyz".index(lx)], B - x0["xyz".index(ly)])
37 u = us[n]
38 far.append(np.abs(u[np.isfinite(u) & (r > np.pi / k)]))
39 finite = np.concatenate(far)
40 umax = np.percentile(finite, 99)
41 ref = np.median(finite)
43 fig, axes = plt.subplots(2, 3, figsize=(13, 8.4), constrained_layout=True)
44 for col, (n, lx, ly, fixed) in enumerate(SLICES):
45 u = us[n]
46 ax = axes[0, col]
47 im0 = ax.imshow(masked(u.real), origin="lower", extent=ext, cmap="RdBu_r",
48 vmin=-umax, vmax=umax, interpolation="bilinear")
49 ax.set_title(f"Re u, {n}-plane at {'xyz'[fixed]} = {x0[fixed]:.2f}")
50 ax = axes[1, col]
51 with np.errstate(divide="ignore", invalid="ignore"):
52 level = 20 * np.log10(np.abs(u) / ref)
53 im1 = ax.imshow(masked(level), origin="lower", extent=ext, cmap="Blues_r",
54 vmin=-30, vmax=15, interpolation="bilinear")
55 ax.set_title(f"level (dB re median), {n}-plane")
56 for ax in axes[:, col]:
57 ax.plot(x0["xyz".index(lx)], x0["xyz".index(ly)], marker="o", ms=5,
58 mfc="none", mec="0.15", mew=1)
59 ax.set_xlabel(lx)
60 ax.set_ylabel(ly)
61 ax.set_aspect("equal")
62 fig.colorbar(im0, ax=axes[0, :], shrink=0.85, label="Re u")
63 fig.colorbar(im1, ax=axes[1, :], shrink=0.85, label="dB")
64 wl = 2 * np.pi / k
65 fig.suptitle(f"Point source in a rounded cubical room: k = {k:g} "
66 f"(room width = {2 * float(d['a']) / wl:.2f} wavelengths), "
67 f"wall admittance {-float(d['lam']):g}")
68 fig.savefig(os.path.join(args.out, "slices.png"), dpi=130)
69 plt.close(fig)
71 # Animation: Re(u e^{-i w t}) over one period on the three planes
72 fig, axes = plt.subplots(1, 3, figsize=(13, 4.6), constrained_layout=True)
73 ims = []
74 for ax, (n, lx, ly, fixed) in zip(axes, SLICES):
75 im = ax.imshow(masked(us[n].real), origin="lower", extent=ext, cmap="RdBu_r",
76 vmin=-umax, vmax=umax, interpolation="bilinear")
77 ax.set_title(f"{n}-plane")
78 ax.set_xlabel(lx)
79 ax.set_ylabel(ly)
80 ims.append(im)
81 fig.colorbar(ims[-1], ax=axes, shrink=0.85, label="Re(u e^{-iωt})")
82 title = fig.suptitle("")
84 def frame(i):
85 ph = np.exp(-2j * np.pi * i / args.frames)
86 for im, (n, *_) in zip(ims, SLICES):
87 im.set_data(masked((us[n] * ph).real))
88 title.set_text(f"pressure over one period, t/T = {i / args.frames:.2f}")
89 return ims
91 anim = FuncAnimation(fig, frame, frames=args.frames, blit=False)
92 anim.save(os.path.join(args.out, "animation.gif"), writer=PillowWriter(fps=12), dpi=80)
93 plt.close(fig)
94 print(f"wrote {args.out}/slices.png and {args.out}/animation.gif")
97if __name__ == "__main__":
98 main()