9b737b8GPU implementation of fmm3dbie's interior Helmholtz impedance solver, with a showcase room problemJeremy Magland 1"""The showcase problem, shared by both backends.
3A point source at x0 inside a rounded cube (the superellipsoid
4|x|^p + |y|^p + |z|^p = a^p) with impedance walls. With time dependence
5exp(-i omega t) the pressure amplitude u satisfies
7 (Delta + k^2) u = -delta(x - x0) in the room,
8 du/dn + i k lambda u = 0 on the walls,
10with n the outward normal and lambda = -beta, where beta > 0 is the wall
11admittance. We write u = G + v with G the free-space field of the source, and
12find v from the right-preconditioned combined field integral equation of
13fmm3dbie (helm_rpcomb_imp_solver, interior mode).
15This module has no dependence on either backend.
16"""
18import json
19import os
20import platform
21import socket
22import subprocess
24import numpy as np
26# Named cases. Both backends discretize the surface identically, so the same
27# case gives the same unknowns on both.
28CASES = {
29 "small": dict(k=10.0, beta=0.2, m=4, norder=6), # 3.2 wavelengths, 4,704 points
30 "medium": dict(k=10.0, beta=0.2, m=6, norder=8), # 3.2 wavelengths, 17,496 points
31 "large": dict(k=40.0, beta=0.05, m=10, norder=8), # 12.7 wavelengths, 48,600 points
32}
33A, P = 1.0, 8 # room half-width and superellipsoid exponent
34SOURCE = np.array([-0.35, 0.2, 0.1]) # the point source
35EXTERIOR_SOURCE = np.array([1.7, 0.6, -0.5]) # for the self-test (outside the room)
36ALPHA = 1.0
39def green(x, x0, k):
40 d = x - x0[:, None]
41 r = np.linalg.norm(d, axis=0)
42 return np.exp(1j * k * r) / (4 * np.pi * r)
45def green_dn(x, nrm, x0, k):
46 d = x - x0[:, None]
47 r = np.linalg.norm(d, axis=0)
48 g = np.exp(1j * k * r) / (4 * np.pi * r)
49 return g * (1j * k - 1 / r) * np.sum(d * nrm, axis=0) / r
52def inside(x, a=A, p=P, margin=0.0):
53 return np.sum(np.abs(x / a) ** p, axis=0) < (1 - margin) ** p
56def room_rhs(xs, nrm, k, lam):
57 """Impedance data for v = u - G: -(dG/dn + i k lambda G)."""
58 return -(green_dn(xs, nrm, SOURCE, k) + 1j * k * lam * green(xs, SOURCE, k))
61def selftest_rhs(xs, nrm, k, lam):
62 """Impedance data of the field of a source outside the room."""
63 return green_dn(xs, nrm, EXTERIOR_SOURCE, k) + 1j * k * lam * green(xs, EXTERIOR_SOURCE, k)
66def selftest_points(n=50):
67 return np.random.default_rng(0).uniform(-0.8 * A, 0.8 * A, size=(3, n))
70def slice_points(ngrid):
71 """Points on three planes through the source; returns grid, {name: (pts, mask, shape)}."""
72 g = np.linspace(-A, A, ngrid)
73 out = {}
74 for name, (i, j, fixed) in {"xy": (0, 1, 2), "xz": (0, 2, 1), "yz": (1, 2, 0)}.items():
75 Ag, Bg = np.meshgrid(g, g, indexing="xy")
76 pts = np.zeros((3, Ag.size))
77 pts[i], pts[j], pts[fixed] = Ag.ravel(), Bg.ravel(), SOURCE[fixed]
78 out[name] = (pts, inside(pts, margin=1e-3), Ag.shape)
79 return g, out
82def energy_balance(k, lam, v_at_source, wts, u_wall):
83 """Power radiated by the source, k/(4 pi) + Im v(x0), and power absorbed by the
84 walls, -k Re(lambda) int |u|^2 dS. They must agree."""
85 p_src = k / (4 * np.pi) + v_at_source.imag
86 p_wall = -k * np.real(lam) * np.sum(wts * np.abs(u_wall) ** 2)
87 return p_src, p_wall
90def machine_info():
91 info = dict(host=socket.gethostname(), platform=platform.platform(), cpu_count=os.cpu_count())
92 try:
93 with open("/proc/cpuinfo") as f:
94 for line in f:
95 if line.startswith("model name"):
96 info["cpu"] = line.split(":", 1)[1].strip()
97 break
98 except OSError:
99 pass
100 try:
101 gpu = subprocess.run(["nvidia-smi", "--query-gpu=name,memory.total", "--format=csv,noheader"],
102 capture_output=True, text=True, timeout=10)
103 if gpu.returncode == 0:
104 info["gpu"] = gpu.stdout.strip().splitlines()[0]
105 except (OSError, subprocess.SubprocessError):
106 pass
107 return info
110def save(outdir, summary, **arrays):
111 os.makedirs(outdir, exist_ok=True)
112 np.savez(os.path.join(outdir, "field.npz"), **arrays)
113 with open(os.path.join(outdir, "summary.json"), "w") as f:
114 json.dump(summary, f, indent=2)