1"""Modal series solution of string-pluck-v1 (see the appendix of the challenge README).
3Usage: python modal_string.py <challenge.json> <out-dir> [N ...]
5Writes modal-n<N>.wav for each truncation order N and prints a JSON list of per-solution
6metadata (timing, parameters) to stdout.
7"""
8from __future__ import annotations
10import json
11import sys
12import time
13from pathlib import Path
15import numpy as np
17sys.path.insert(0, str(Path(__file__).resolve().parents[1]))
18from common import hardware_string, load_challenge, write_wav # noqa: E402
19from fields import field_grid, write_fields # noqa: E402
22def modal_coefficients(p: dict, n_modes: int):
23 n = np.arange(1, n_modes + 1)
24 k = n * np.pi / p["L"]
25 q = np.pi / p["w"]
26 A = (2.0 * p["h"] / p["L"]) * np.sin(k * p["x_p"]) * np.sin(k * p["w"]) * q**2 / (k * (q**2 - k**2))
27 alpha = p["sigma0"] + 0.5 * p["kappa"] * k**2
28 omega2 = (p["c"] * k) ** 2 - alpha**2
29 if np.any(omega2 <= 0):
30 raise ValueError("an overdamped mode is within the requested truncation; not handled")
31 omega = np.sqrt(omega2)
32 return k, A, alpha, omega
35def pickup_velocity(p: dict, n_modes: int, t: np.ndarray) -> np.ndarray:
36 k, A, alpha, omega = modal_coefficients(p, n_modes)
37 coeff = -A * (p["c"] * k) ** 2 / omega * np.sin(k * p["x_m"]) # per-mode amplitude of a_n'(t) sin(k_n x_m)
38 s = np.zeros_like(t)
39 block = 64 # modes per block, to bound memory for large N
40 for i in range(0, n_modes, block):
41 sl = slice(i, i + block)
42 e = np.exp(-np.outer(t, alpha[sl]))
43 s += (e * np.sin(np.outer(t, omega[sl]))) @ coeff[sl]
44 return s
47def displacement_field(p: dict, n_modes: int, x: np.ndarray, t: np.ndarray) -> np.ndarray:
48 """u(x_i, t_j) as an array of shape (nt, nx)."""
49 k, A, alpha, omega = modal_coefficients(p, n_modes)
50 u = np.zeros((len(t), len(x)))
51 block = 64
52 for i in range(0, n_modes, block):
53 sl = slice(i, i + block)
54 e = np.exp(-np.outer(t, alpha[sl]))
55 a = A[sl] * e * (np.cos(np.outer(t, omega[sl])) + (alpha[sl] / omega[sl]) * np.sin(np.outer(t, omega[sl])))
56 u += a @ np.sin(np.outer(k[sl], x))
57 return u
60def main():
61 challenge = load_challenge(Path(sys.argv[1]))
62 out_dir = Path(sys.argv[2])
63 args = sys.argv[3:]
64 fields_for = set()
65 if "--fields" in args:
66 i = args.index("--fields")
67 fields_for = {int(v) for v in args[i + 1].split(",")}
68 args = args[:i] + args[i + 2:]
69 orders = [int(a) for a in args] or [4, 8, 16, 32, 64, 256, 1024]
70 p = challenge["parameters"]
71 out = challenge["output"]
72 fs = out["sample_rate_hz"]
73 t = np.arange(out["num_samples"]) / fs
74 meta = []
75 for n in orders:
76 pickup_velocity(p, min(n, 8), t[:16]) # warm up
77 t0 = time.perf_counter()
78 s = pickup_velocity(p, n, t)
79 dt = time.perf_counter() - t0
80 write_wav(out_dir / f"modal-n{n}.wav", s, challenge)
81 entry = {
82 "solution_id": f"modal-n{n}",
83 "output": f"{out_dir.name}/modal-n{n}.wav",
84 "method": "modal series, truncated",
85 "parameters": {"n_modes": n},
86 "timing_seconds": round(dt, 4),
87 "hardware": hardware_string(),
88 "peak_si": float(np.max(np.abs(s))),
89 }
90 if n in fields_for:
91 x, tt = field_grid(challenge)
92 write_fields(out_dir / f"modal-n{n}.h5", x, tt, displacement_field(p, n, x, tt), challenge, f"modal-n{n}")
93 entry["fields"] = {"path": f"{out_dir.name}/modal-n{n}.h5", "format": "hdf5"}
94 meta.append(entry)
95 print(f"N={n:5d} peak={np.max(np.abs(s)):.4f} m/s t={dt:.3f}s", file=sys.stderr)
96 print(json.dumps(meta, indent=1))
99if __name__ == "__main__":
100 main()