/ challenges-old / string-solutions
Sign in
challenges-old / string-solutions
main
Branches
main
string-solutions / solvers / fd / fd_string.py
134 lines · 5.5 KBCodeBlameHistory
1c822d5Add modal-series and fd-second-order submissions for string-pluck-v1, with results workflowJeremy Magland 1"""Finite-difference solution of string-pluck-v1, second order in space and time.
3With D the standard three-point second difference (zero at the fixed ends) and u^n the grid
4values at step n, the scheme is
6 (u^{n+1} - 2u^n + u^{n-1})/dt^2 = c^2 D u^n / dx^2
7 - sigma0 (u^{n+1} - u^{n-1}) / dt
8 + kappa D (u^{n+1} - u^{n-1}) / (2 dt dx^2).
10The Kelvin-Voigt term is centered in time, so the scheme is implicit in the damping and a
11tridiagonal system is solved at every step. (An earlier explicit variant with a backward
12difference for that term is only first order in time; its error cancelled against the
13second-order dispersion error near m = 8 and the sequence was not monotone.) Stability
14requires the Courant number c dt/dx <= 1.
16The time step is dt = 1/(fs*m) for an integer oversampling factor m, and the grid spacing is
17dx = c dt / lam with lam = 0.9, so the scheme is deliberately dispersive. The pickup velocity is
18(u^{n+1} - u^{n-1})/(2 dt) linearly interpolated to x_m, taken every m-th step.
20Usage: python fd_string.py <challenge.json> <out-dir> [m ...]
21"""
22from __future__ import annotations
24import json
25import sys
26import time
27from pathlib import Path
29import numpy as np
30from scipy.linalg import solve_banded
32sys.path.insert(0, str(Path(__file__).resolve().parents[1]))
33from common import hardware_string, initial_shape, load_challenge, write_wav # noqa: E402
36LAM = 0.9
8f29519Export displacement fields for the chain tops, add this repository's own site with a field viewer, and leave the comparison to the challenge repositoryJeremy Magland 39def solve(p: dict, fs: int, num_samples: int, m: int, lam: float = LAM, field_x=None, field_t=None):
40 """Returns (pickup samples, info, fields) where fields is u on (field_t, field_x) or None.
42 Field frames are taken at the step nearest each requested time and linearly interpolated in x."""
44 n_int = int(np.floor(p["L"] * lam / (p["c"] * dt))) # intervals; floor keeps dx >= c dt / lam
45 dx = p["L"] / n_int
46 x = np.linspace(0.0, p["L"], n_int + 1)
47 lam2 = (p["c"] * dt / dx) ** 2
48 s = p["sigma0"] * dt
49 beta = p["kappa"] * dt / (2.0 * dx**2)
51 def D(u):
52 out = np.zeros_like(u)
53 out[1:-1] = u[:-2] - 2.0 * u[1:-1] + u[2:]
54 return out
56 # Left-hand operator on interior points: (1 + s) I - beta D, tridiagonal, in solve_banded form.
57 n_in = n_int - 1
58 ab = np.zeros((3, n_in))
59 ab[0, 1:] = -beta
60 ab[1, :] = (1.0 + s) + 2.0 * beta
61 ab[2, :-1] = -beta
63 u_prev = initial_shape(x, p)
64 u_cur = u_prev + 0.5 * dt**2 * (p["c"] ** 2) * D(u_prev) / dx**2 # zero initial velocity, second-order start
65 xi = p["x_m"] / dx
66 i0 = int(np.floor(xi))
67 w1 = xi - i0
68 out = np.zeros(num_samples)
69 total_steps = num_samples * m
70 u_next = np.zeros_like(u_cur)
72 frame_steps = {}
73 if field_x is not None:
74 fields = np.zeros((len(field_t), len(field_x)))
75 for j, tj in enumerate(field_t):
76 frame_steps.setdefault(int(round(tj / dt)), []).append(j)
77 for j in frame_steps.get(0, []):
78 fields[j] = np.interp(field_x, x, u_prev)
79 for j in frame_steps.get(1, []):
80 fields[j] = np.interp(field_x, x, u_cur)
82 rhs = 2.0 * u_cur + lam2 * D(u_cur) - (1.0 - s) * u_prev - beta * D(u_prev)
83 u_next[1:-1] = solve_banded((1, 1), ab, rhs[1:-1], overwrite_b=True, check_finite=False)
84 if n % m == 0:
85 k = n // m
86 if k < num_samples:
87 v = (u_next - u_prev) / (2.0 * dt)
88 out[k] = (1.0 - w1) * v[i0] + w1 * v[i0 + 1]
90 fields[j] = np.interp(field_x, x, u_next)
1c822d5Add modal-series and fd-second-order submissions for string-pluck-v1, with results workflowJeremy Magland 91 u_prev, u_cur, u_next = u_cur, u_next, u_prev
8f29519Export displacement fields for the chain tops, add this repository's own site with a field viewer, and leave the comparison to the challenge repositoryJeremy Magland 92 return out, {"m": m, "dt": dt, "dx": dx, "n_points": n_int + 1, "courant": p["c"] * dt / dx}, fields
95def main():
96 challenge = load_challenge(Path(sys.argv[1]))
97 out_dir = Path(sys.argv[2])
99 fields_for = set()
100 if "--fields" in args:
101 i = args.index("--fields")
102 fields_for = {int(v) for v in args[i + 1].split(",")}
103 args = args[:i] + args[i + 2:]
104 ms = [int(a) for a in args] or [1, 2, 4, 8, 16]
106 out = challenge["output"]
107 meta = []
108 for m in ms:
109 t0 = time.perf_counter()
112 sid = f"fd-m{m}"
113 write_wav(out_dir / f"{sid}.wav", sig, challenge)
116 "output": f"{out_dir.name}/{sid}.wav",
117 "method": "finite differences, second order, implicit damping",
118 "parameters": {"oversampling": m, "n_points": info["n_points"], "dx_m": info["dx"], "dt_s": info["dt"], "courant": round(info["courant"], 6)},
119 "timing_seconds": round(dt, 3),
120 "hardware": hardware_string(),
121 "peak_si": float(np.max(np.abs(sig))),
123 if m in fields_for:
124 x, tt = field_grid(challenge)
125 _, _, fld = solve(p, out["sample_rate_hz"], out["num_samples"], m, field_x=x, field_t=tt)
126 write_fields(out_dir / f"{sid}.h5", x, tt, fld, challenge, sid)
127 entry["fields"] = {"path": f"{out_dir.name}/{sid}.h5", "format": "hdf5"}
128 meta.append(entry)
1c822d5Add modal-series and fd-second-order submissions for string-pluck-v1, with results workflowJeremy Magland 129 print(f"m={m} points={info['n_points']} peak={np.max(np.abs(sig)):.4f} m/s t={dt:.2f}s", file=sys.stderr)
130 print(json.dumps(meta, indent=1))
133if __name__ == "__main__":
134 main()
↑↓ move↵openescclose