/ proof-of-concept / fmm3dbie-gpu-poc
Sign in
proof-of-concept / fmm3dbie-gpu-poc
bb9c969
Branches
main
fmm3dbie-gpu-poc / plot.py
191 lines · 8.4 KBCodeBlameHistory
bb9c969Add the cymbal problem: a listener near a rigid cymbal in free spaceJeremy Magland 1"""Plot a 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. Works for the room
3and the cymbal problems."""
5import argparse
9import matplotlib
11matplotlib.use("Agg")
12import matplotlib.pyplot as plt
13import numpy as np
14from matplotlib.animation import FuncAnimation, PillowWriter
16SLICES = [("xy", "x", "y", 2), ("xz", "x", "z", 1), ("yz", "y", "z", 0)]
19def masked(a):
20 return np.ma.masked_invalid(a)
23def main():
24 ap = argparse.ArgumentParser()
25 ap.add_argument("out", help="a results directory written by run.py")
26 ap.add_argument("--frames", type=int, default=32)
27 args = ap.parse_args()
29 d = np.load(os.path.join(args.out, "field.npz"))
bb9c969Add the cymbal problem: a listener near a rigid cymbal in free spaceJeremy Magland 30 with open(os.path.join(args.out, "summary.json")) as f:
31 if json.load(f).get("problem", "room") == "cymbal":
32 return plot_cymbal(d, args)
34 ext = [g[0], g[-1], g[0], g[-1]]
35 us = {n: d[f"u_{n}"] for n, *_ in SLICES}
37 # Color limits from the field more than half a wavelength from the source
38 A, B = np.meshgrid(g, g, indexing="xy")
39 far = []
40 for n, lx, ly, fixed in SLICES:
41 r = np.hypot(A - x0["xyz".index(lx)], B - x0["xyz".index(ly)])
42 u = us[n]
43 far.append(np.abs(u[np.isfinite(u) & (r > np.pi / k)]))
44 finite = np.concatenate(far)
45 umax = np.percentile(finite, 99)
46 ref = np.median(finite)
48 fig, axes = plt.subplots(2, 3, figsize=(13, 8.4), constrained_layout=True)
49 for col, (n, lx, ly, fixed) in enumerate(SLICES):
50 u = us[n]
51 ax = axes[0, col]
52 im0 = ax.imshow(masked(u.real), origin="lower", extent=ext, cmap="RdBu_r",
53 vmin=-umax, vmax=umax, interpolation="bilinear")
54 ax.set_title(f"Re u, {n}-plane at {'xyz'[fixed]} = {x0[fixed]:.2f}")
55 ax = axes[1, col]
56 with np.errstate(divide="ignore", invalid="ignore"):
57 level = 20 * np.log10(np.abs(u) / ref)
58 im1 = ax.imshow(masked(level), origin="lower", extent=ext, cmap="Blues_r",
59 vmin=-30, vmax=15, interpolation="bilinear")
60 ax.set_title(f"level (dB re median), {n}-plane")
61 for ax in axes[:, col]:
62 ax.plot(x0["xyz".index(lx)], x0["xyz".index(ly)], marker="o", ms=5,
63 mfc="none", mec="0.15", mew=1)
64 ax.set_xlabel(lx)
65 ax.set_ylabel(ly)
66 ax.set_aspect("equal")
67 fig.colorbar(im0, ax=axes[0, :], shrink=0.85, label="Re u")
68 fig.colorbar(im1, ax=axes[1, :], shrink=0.85, label="dB")
69 wl = 2 * np.pi / k
70 fig.suptitle(f"Point source in a rounded cubical room: k = {k:g} "
71 f"(room width = {2 * float(d['a']) / wl:.2f} wavelengths), "
72 f"wall admittance {-float(d['lam']):g}")
73 fig.savefig(os.path.join(args.out, "slices.png"), dpi=130)
74 plt.close(fig)
76 # Animation: Re(u e^{-i w t}) over one period on the three planes
77 fig, axes = plt.subplots(1, 3, figsize=(13, 4.6), constrained_layout=True)
78 ims = []
79 for ax, (n, lx, ly, fixed) in zip(axes, SLICES):
80 im = ax.imshow(masked(us[n].real), origin="lower", extent=ext, cmap="RdBu_r",
81 vmin=-umax, vmax=umax, interpolation="bilinear")
82 ax.set_title(f"{n}-plane")
83 ax.set_xlabel(lx)
84 ax.set_ylabel(ly)
85 ims.append(im)
86 fig.colorbar(ims[-1], ax=axes, shrink=0.85, label="Re(u e^{-iωt})")
87 title = fig.suptitle("")
89 def frame(i):
90 ph = np.exp(-2j * np.pi * i / args.frames)
91 for im, (n, *_) in zip(ims, SLICES):
92 im.set_data(masked((us[n] * ph).real))
93 title.set_text(f"pressure over one period, t/T = {i / args.frames:.2f}")
94 return ims
96 anim = FuncAnimation(fig, frame, frames=args.frames, blit=False)
97 anim.save(os.path.join(args.out, "animation.gif"), writer=PillowWriter(fps=12), dpi=80)
98 plt.close(fig)
99 print(f"wrote {args.out}/slices.png and {args.out}/animation.gif")
bb9c969Add the cymbal problem: a listener near a rigid cymbal in free spaceJeremy Magland 102CYMBAL_SLICES = ["xz", "yz", "xy"]
105def plot_cymbal(d, args):
106 """Slices through the evaluation box around the cymbal. u(x) is the pressure at the
107 listener due to a unit source at x (reciprocity). The level is shown relative to
108 the free field, 20 log10 |u / G|, which is 0 dB where the cymbal has no effect."""
109 x0, k, freq = d["x0"], float(d["k"]), float(d["freq"])
110 rz, R = d["outline"], float(d["radius"])
111 names = [n for n in CYMBAL_SLICES if f"u_{n}" in d]
112 panels = []
113 for n in names:
114 ga, gb = d[f"axes_{n}"]
115 fixed, val = d[f"plane_{n}"]
116 fixed = int(fixed)
117 ia, ib = [i for i in range(3) if i != fixed]
118 A, B = np.meshgrid(ga, gb, indexing="xy")
119 pts = np.zeros((3, A.size))
120 pts[ia], pts[ib], pts[fixed] = A.ravel(), B.ravel(), val
121 r = np.linalg.norm(pts - x0[:, None], axis=0)
122 G = (np.exp(1j * k * r) / (4 * np.pi * r)).reshape(A.shape)
123 panels.append((n, ga, gb, ia, ib, fixed, val, d[f"u_{n}"], G, r.reshape(A.shape)))
124 far = np.concatenate([np.abs(u[np.isfinite(u) & (r > np.pi / k)]) for *_, u, G, r in panels])
125 umax = np.percentile(far, 99)
127 def outline(ax, ia, ib, fixed, val):
128 if fixed == 2: # horizontal plane: the cymbal's rim, projected
129 t = np.linspace(0, 2 * np.pi, 200)
130 ax.plot(R * np.cos(t), R * np.sin(t), color="0.2", lw=0.8, ls="--")
131 elif abs(val) < R: # vertical plane: the cross section, at distance |val| from the axis
132 r, z = rz
133 m = r >= abs(val)
134 s = np.sqrt(np.maximum(r[m] ** 2 - val ** 2, 0))
135 for sign in (1, -1):
136 ax.fill(sign * s, z[m], color="0.25", lw=0)
138 fig, axes = plt.subplots(2, len(panels), figsize=(4.4 * len(panels), 8.4), constrained_layout=True,
139 squeeze=False)
140 for col, (n, ga, gb, ia, ib, fixed, val, u, G, r) in enumerate(panels):
141 ext = [ga[0], ga[-1], gb[0], gb[-1]]
142 ax = axes[0, col]
143 im0 = ax.imshow(masked(u.real), origin="lower", extent=ext, cmap="RdBu_r",
144 vmin=-umax, vmax=umax, interpolation="bilinear")
145 ax.set_title(f"Re u, {n}-plane at {'xyz'[fixed]} = {val:.2f} m")
146 ax = axes[1, col]
147 with np.errstate(divide="ignore", invalid="ignore"):
148 level = 20 * np.log10(np.abs(u) / np.abs(G))
149 im1 = ax.imshow(masked(level), origin="lower", extent=ext, cmap="PuOr_r",
150 vmin=-15, vmax=15, interpolation="bilinear")
151 ax.set_title(f"level re free field (dB), {n}-plane")
152 for ax in axes[:, col]:
153 outline(ax, ia, ib, fixed, val)
154 if abs(x0[fixed] - val) < 1e-9:
155 ax.plot(x0[ia], x0[ib], marker="o", ms=5, mfc="none", mec="0.15", mew=1)
156 ax.set_xlabel(f"{'xyz'[ia]} (m)")
157 ax.set_ylabel(f"{'xyz'[ib]} (m)")
158 ax.set_aspect("equal")
159 fig.colorbar(im0, ax=axes[0, :], shrink=0.85, label="Re u")
160 fig.colorbar(im1, ax=axes[1, :], shrink=0.85, label="dB")
161 wl = 2 * np.pi / k
162 fig.suptitle(f"Listener at ({x0[0]:.2f}, {x0[1]:.2f}, {x0[2]:.2f}) m near a rigid cymbal, "
163 f"{freq:g} Hz (diameter = {2 * R / wl:.2f} wavelengths)\nBy reciprocity u(x) is the "
164 f"pressure at the listener due to a unit source at x")
165 fig.savefig(os.path.join(args.out, "slices.png"), dpi=130)
166 plt.close(fig)
168 # Animation of the vertical plane through the listener
169 n, ga, gb, ia, ib, fixed, val, u, G, r = panels[0]
170 fig, ax = plt.subplots(figsize=(6, 5.4), constrained_layout=True)
171 im = ax.imshow(masked(u.real), origin="lower", extent=[ga[0], ga[-1], gb[0], gb[-1]], cmap="RdBu_r",
172 vmin=-umax, vmax=umax, interpolation="bilinear")
173 outline(ax, ia, ib, fixed, val)
174 ax.set_xlabel(f"{'xyz'[ia]} (m)")
175 ax.set_ylabel(f"{'xyz'[ib]} (m)")
176 fig.colorbar(im, ax=ax, shrink=0.85, label="Re(u e^{-iωt})")
177 title = ax.set_title("")
179 def frame(i):
180 im.set_data(masked((u * np.exp(-2j * np.pi * i / args.frames)).real))
181 title.set_text(f"{n}-plane, t/T = {i / args.frames:.2f}")
182 return [im]
184 anim = FuncAnimation(fig, frame, frames=args.frames, blit=False)
185 anim.save(os.path.join(args.out, "animation.gif"), writer=PillowWriter(fps=12), dpi=80)
186 plt.close(fig)
187 print(f"wrote {args.out}/slices.png and {args.out}/animation.gif")
191 main()
↑↓ move↵openescclose