bb9c969Add the cymbal problem: a listener near a rigid cymbal in free spaceJeremy Magland 1"""Helpers shared by the problems: the Green's function, boundary data, machine
2information and saving results. No dependence on either backend."""
4import json
5import os
6import platform
7import socket
8import subprocess
10import numpy as np
13def green(x, x0, k):
14 d = x - x0[:, None]
15 r = np.linalg.norm(d, axis=0)
16 return np.exp(1j * k * r) / (4 * np.pi * r)
19def green_dn(x, nrm, x0, k):
20 d = x - x0[:, None]
21 r = np.linalg.norm(d, axis=0)
22 g = np.exp(1j * k * r) / (4 * np.pi * r)
23 return g * (1j * k - 1 / r) * np.sum(d * nrm, axis=0) / r
26def impedance_data(xs, nrm, x0, k, lam):
27 """du/dn + i k lambda u on the surface for u the field of a unit source at x0."""
28 return green_dn(xs, nrm, x0, k) + 1j * k * lam * green(xs, x0, k)
31def source_sum(x, x0s, ws, k, nrm=None, lam=0.0, chunk=64, xp=np):
32 """Field sum_j ws[j] G(x, x0s[:, j]) at points x (3, n), or with nrm given its
33 impedance data du/dn + i k lambda u. xp is numpy or cupy; the result is NumPy."""
34 x, x0s, ws = xp.asarray(x), xp.asarray(x0s), xp.asarray(ws)
35 nrm = None if nrm is None else xp.asarray(nrm)
36 out = xp.zeros(x.shape[1], dtype=complex)
37 for i in range(0, x0s.shape[1], chunk):
38 d = x[:, None, :] - x0s[:, i:i + chunk, None] # (3, chunk, n)
39 r = xp.sqrt(xp.sum(d * d, axis=0))
40 g = xp.exp(1j * k * r) / (4 * np.pi * r)
41 if nrm is not None:
42 g = g * ((1j * k - 1 / r) * xp.einsum("icn,in->cn", d, nrm) / r + 1j * k * lam)
43 out += ws[i:i + chunk] @ g
44 return out.get() if xp is not np else out
47def machine_info():
48 info = dict(host=socket.gethostname(), platform=platform.platform(), cpu_count=os.cpu_count())
49 try:
50 with open("/proc/cpuinfo") as f:
51 for line in f:
52 if line.startswith("model name"):
53 info["cpu"] = line.split(":", 1)[1].strip()
54 break
55 except OSError:
56 pass
57 try:
58 gpu = subprocess.run(["nvidia-smi", "--query-gpu=name,memory.total", "--format=csv,noheader"],
59 capture_output=True, text=True, timeout=10)
60 if gpu.returncode == 0:
61 info["gpu"] = gpu.stdout.strip().splitlines()[0]
62 except (OSError, subprocess.SubprocessError):
63 pass
64 return info
67def save(outdir, summary, **arrays):
68 os.makedirs(outdir, exist_ok=True)
69 np.savez(os.path.join(outdir, "field.npz"), **arrays)
70 with open(os.path.join(outdir, "summary.json"), "w") as f:
71 json.dump(summary, f, indent=2)