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, pluck_load_q, 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 Q = pluck_load_q(p)
27 A = (2.0 * Q / (p["L"] * p["w"] * k**3)) * np.sin(k * p["x_p"]) * np.sin(k * p["w"]) * q**2 / (q**2 - k**2)
28 alpha = p["sigma0"] + 0.5 * p["kappa"] * k**2
29 omega2 = (p["c"] * k) ** 2 - alpha**2
30 if np.any(omega2 <= 0):
31 raise ValueError("an overdamped mode is within the requested truncation; not handled")
32 omega = np.sqrt(omega2)
33 return k, A, alpha, omega
36def pickup_velocity(p: dict, n_modes: int, t: np.ndarray) -> np.ndarray:
37 k, A, alpha, omega = modal_coefficients(p, n_modes)
38 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)
39 s = np.zeros_like(t)
40 block = 64 # modes per block, to bound memory for large N
41 for i in range(0, n_modes, block):
42 sl = slice(i, i + block)
43 e = np.exp(-np.outer(t, alpha[sl]))
44 s += (e * np.sin(np.outer(t, omega[sl]))) @ coeff[sl]
45 return s
48def displacement_field(p: dict, n_modes: int, x: np.ndarray, t: np.ndarray) -> np.ndarray:
49 """u(x_i, t_j) as an array of shape (nt, nx)."""
50 k, A, alpha, omega = modal_coefficients(p, n_modes)
51 u = np.zeros((len(t), len(x)))
52 block = 64
53 for i in range(0, n_modes, block):
54 sl = slice(i, i + block)
55 e = np.exp(-np.outer(t, alpha[sl]))
56 a = A[sl] * e * (np.cos(np.outer(t, omega[sl])) + (alpha[sl] / omega[sl]) * np.sin(np.outer(t, omega[sl])))
57 u += a @ np.sin(np.outer(k[sl], x))
58 return u
61def main():
62 challenge = load_challenge(Path(sys.argv[1]))
63 out_dir = Path(sys.argv[2])
64 args = sys.argv[3:]
65 fields_for = set()
66 if "--fields" in args:
67 i = args.index("--fields")
68 fields_for = {int(v) for v in args[i + 1].split(",")}
69 args = args[:i] + args[i + 2:]
70 orders = [int(a) for a in args] or [4, 8, 16, 32, 64, 256, 1024]
71 p = challenge["parameters"]
72 out = challenge["output"]
73 fs = out["sample_rate_hz"]
74 t = np.arange(out["num_samples"]) / fs
75 meta = []
76 for n in orders:
77 pickup_velocity(p, min(n, 8), t[:16]) # warm up
78 t0 = time.perf_counter()
79 s = pickup_velocity(p, n, t)
80 dt = time.perf_counter() - t0
81 write_wav(out_dir / f"modal-n{n}.wav", s, challenge)
82 entry = {
83 "solution_id": f"modal-n{n}",
84 "output": f"{out_dir.name}/modal-n{n}.wav",
85 "method": "modal series, truncated",
86 "parameters": {"n_modes": n},
87 "timing_seconds": round(dt, 4),
88 "hardware": hardware_string(),
89 "peak_si": float(np.max(np.abs(s))),
90 }
91 if n in fields_for:
92 x, tt = field_grid(challenge)
93 write_fields(out_dir / f"modal-n{n}.h5", x, tt, displacement_field(p, n, x, tt), challenge, f"modal-n{n}")
94 entry["fields"] = {"path": f"{out_dir.name}/modal-n{n}.h5", "format": "hdf5"}
95 meta.append(entry)
96 print(f"N={n:5d} peak={np.max(np.abs(s)):.4f} m/s t={dt:.3f}s", file=sys.stderr)
97 print(json.dumps(meta, indent=1))
100if __name__ == "__main__":
101 main()