Strike schedules for experiments
3 changed files+53−10
solvers/vonkarman-modal-disk/__pycache__/bfs_disk.cpython-313.pycdeleted+0−0View file
Binary file not shown.
solvers/vonkarman-modal-disk/__pycache__/vk_modal_disk.cpython-313.pycdeleted+0−0View file
Binary file not shown.
solvers/vonkarman-modal-disk/vk_modal_disk.pymodified+53−10View file
@@ -44,7 +44,9 @@ neglects the mass and stiffness of a 12 mm disk at the apex; a stated approximat
4444 With --device cuda the coupling tensor is assembled and the time stepping runs on a GPU through PyTorch
4545 (the modes and quadrature stay in numpy/scipy); the CPU path is the reference and the two agree to round-off.
4646
47-Usage: python vk_modal_disk.py <challenge.json> <out-dir> [--fields tag,...] [--force F] [--duration T] [--seconds S] [--full] [--device cuda] [tag ...]
47+Usage: python vk_modal_disk.py <challenge.json> <out-dir> [--fields tag,...] [--force F] [--duration T] [--seconds S] [--full] [--device cuda] [--strikes FILE] [tag ...]
48+ --strikes FILE replaces the challenge's strike by a schedule, a JSON list of {"t", "force", "duration", "x", "y"}, for
49+ experiments outside the challenge (the output is then not a valid submission).
4850 where a tag is g<grid spacing in mm>-w<transverse modes>-a<airy modes>-m<oversampling>, e.g. g8-w200-a200-m1
4951 """
5052 from __future__ import annotations
@@ -174,8 +176,11 @@ Paraboloid = Profile # former name
174176
175177
176178 class VonKarmanDisk:
177- def __init__(self, ch: dict, spacing: float, n_w: int, n_airy: int, force_override: float | None = None, log=None, symmetric_y: bool | None = None, device: str = "cpu", duration_override: float | None = None):
179+ def __init__(self, ch: dict, spacing: float, n_w: int, n_airy: int, force_override: float | None = None, log=None, symmetric_y: bool | None = None, device: str = "cpu", duration_override: float | None = None, strikes=None):
180+ """strikes: optional list of {"t", "force", "duration", "x", "y"} replacing the challenge's single strike at
181+ t = 0 (for experiments outside the challenge; the output is then not a valid submission)."""
178182 self.ch = ch
183+ self.strikes = strikes
179184 self.device = device
180185 self.torch = None
181186 if device != "cpu":
@@ -258,6 +263,38 @@ class VonKarmanDisk:
258263 # Rayleigh integral: quadrature weights, delays, modal shapes at the points
259264 _, _, self.tau, self.cq = pl.listener_quadrature(3)
260265 self.PhiQ = Bn @ V # (Q, n_w)
266+ # the strike schedule: (t0, duration, modal force vector) per hit; default is the challenge's single strike at t = 0
267+ if self.strikes:
268+ self.schedule = [(float(st["t"]), float(st["duration"]), self.strike_force(float(st["x"]), float(st["y"]), float(st["force"]))) for st in self.strikes]
269+ else:
270+ self.schedule = [(0.0, p["strike_duration"], self.f)]
271+
272+ def strike_force(self, xs: float, ys: float, F: float) -> np.ndarray:
273+ """Modal force vector of a raised-cosine strike of peak total force F centered at (xs, ys) (horizontal
274+ coordinates of the top face), for the flat disk or the shell, normalized by the surface integral of the bump."""
275+ from numpy.polynomial import legendre as leg
276+
277+ p, pl, pr = self.p, self.plate, self.profile
278+ Rs = p["strike_radius"]
279+ gr, gw = leg.leggauss(32)
280+ rq = (Rs * 1.1) * (gr + 1.0) / 2.0
281+ wr = gw * (Rs * 1.1) / 2.0
282+ nth = 96
283+ th = 2 * np.pi * np.arange(nth) / nth
284+ rr, tt = np.meshgrid(rq, th, indexing="ij")
285+ wa = (wr[:, None] * rr * (2 * np.pi / nth)).ravel()
286+ px = (xs + rr * np.cos(tt)).ravel()
287+ py = (ys + rr * np.sin(tt)).ravel()
288+ if pr is None:
289+ tx, ty, nz = px, py, np.ones_like(px)
290+ else:
291+ tx, ty, _ = pr.top(px, py)
292+ _, _, nz = pr.normal(px, py)
293+ rprime = np.sqrt((tx - xs) ** 2 + (ty - ys) ** 2)
294+ phi = np.where(rprime <= Rs, 0.5 * (1.0 + np.cos(np.pi * np.minimum(rprime / Rs, 1.0))), 0.0)
295+ C = np.sum(wa * phi / nz) # surface integral of the bump over the top face
296+ Phi_s = pl.basis_at(px, py) @ self.V
297+ return -(F / C) * ((wa * phi / nz) @ Phi_s)
261298
262299 def _apply_profile(self, xq, yq, WQ, PsiQ, Pxx, Pyy, Pxy):
263300 """Marguerre shallow-shell terms for the initial shape: the G tensor, the shell modes (rotation T of the
@@ -390,8 +427,10 @@ class VonKarmanDisk:
390427 report = max(1, nsteps // 10)
391428 for nstep in range(nsteps + 1):
392429 t = nstep * dt
393- psi = 0.5 * (1.0 - np.cos(Om * t)) if 0.0 <= t <= T else 0.0
394- Fext = self.nonlinear_force(q) + f * psi
430+ Fext = self.nonlinear_force(q)
431+ for t0, Ts, fvec in self.schedule:
432+ if t0 <= t <= t0 + Ts:
433+ Fext = Fext + fvec * (0.5 * (1.0 - np.cos(2 * np.pi * (t - t0) / Ts)))
395434 q_next = a1 * q + a2 * q_prev + b * Fext
396435 qdd = Fext - w2 * q - 2.0 * al * (q_next - q_prev) / (2.0 * dt)
397436 ring[nstep % (dmax + 1)] = qdd
@@ -421,7 +460,8 @@ class VonKarmanDisk:
421460 dev = t.device(self.device)
422461 T64 = lambda a: t.as_tensor(np.asarray(a, dtype=np.float64), device=dev)
423462 dt, T, Om, nsteps, dmax, fld = pr["dt"], pr["T"], pr["Om"], pr["nsteps"], pr["dmax"], pr["fld"]
424- a1, a2, b, w2, al, G, f, coef = (T64(x) for x in (pr["a1"], pr["a2"], pr["b"], pr["w2"], pr["al"], pr["G"], self.f, self.coef))
463+ a1, a2, b, w2, al, G, coef = (T64(x) for x in (pr["a1"], pr["a2"], pr["b"], pr["w2"], pr["al"], pr["G"], self.coef))
464+ schedule = [(t0, Ts, T64(fvec)) for t0, Ts, fvec in self.schedule]
425465 Hf = self.Hflat_t
426466 Gt = T64(self.G) if self.G is not None else None
427467 L = dmax + 1
@@ -437,15 +477,15 @@ class VonKarmanDisk:
437477 report = max(1, nsteps // 10)
438478 for nstep in range(nsteps + 1):
439479 tt = nstep * dt
440- psi = 0.5 * (1.0 - np.cos(Om * tt)) if 0.0 <= tt <= T else 0.0
441480 v = (Hf @ q).reshape(self.n_airy, self.n_w)
442481 sq = v @ q
443482 if Gt is None:
444483 Fext = v.T @ (coef * sq)
445484 else:
446485 Fext = v.T @ (coef * (sq + 2.0 * (Gt @ q))) + Gt.T @ (coef * sq)
447- if psi:
448- Fext = Fext + f * psi
486+ for t0, Ts, fvec in schedule:
487+ if t0 <= tt <= t0 + Ts:
488+ Fext = Fext + fvec * (0.5 * (1.0 - np.cos(2 * np.pi * (tt - t0) / Ts)))
449489 q_next = a1 * q + a2 * q_prev + b * Fext
450490 qdd = Fext - w2 * q - al * ((q_next - q_prev) * invdt)
451491 ring[nstep % L] = qdd
@@ -517,6 +557,8 @@ def main():
517557 fields_for = v or set()
518558 force = take("--force", float)
519559 duration = take("--duration", float)
560+ strikes_file = take("--strikes", str)
561+ strikes = json.load(open(strikes_file)) if strikes_file else None
520562 seconds = take("--seconds", float)
521563 device = take("--device", str) or "cpu"
522564 full = "--full" in args
@@ -532,7 +574,7 @@ def main():
532574 spacing, n_w, n_a, m = parse_tag(tag)
533575 log(f"{tag}: force={force if force is not None else ch['parameters']['strike_peak_force']} N, duration={duration if duration is not None else ch['parameters']['strike_duration']} s")
534576 t0 = time.perf_counter()
535- model = VonKarmanDisk(ch, spacing, n_w, n_a, force, log=log, symmetric_y=False if full else None, device=device, duration_override=duration)
577+ model = VonKarmanDisk(ch, spacing, n_w, n_a, force, log=log, symmetric_y=False if full else None, device=device, duration_override=duration, strikes=strikes)
536578 t1 = time.perf_counter()
537579 grid = None
538580 if tag in fields_for:
@@ -542,7 +584,7 @@ def main():
542584 grid = (np.linspace(-R, R, nx), np.linspace(-R, R, ny), np.arange(nt) / fg["frame_rate_hz"])
543585 s, fields = model.run(fs, nsamp, m, grid, ch["fields"]["frame_rate_hz"] if grid else None)
544586 dt_total = time.perf_counter() - t0
545- sid = f"vk-{tag}"
587+ sid = f"vk-{tag}" + ("-" + Path(strikes_file).stem if strikes_file else "")
546588 out_dir.mkdir(parents=True, exist_ok=True)
547589 scale = out["signal_scale"] if out.get("signal_scale") else 1.0 # a draft challenge may have no scale yet
548590 wavfile.write(out_dir / f"{sid}.wav", fs, (s / scale).astype(np.float32))
@@ -566,6 +608,7 @@ def main():
566608 "hardware": hardware_string(device),
567609 "peak_si": float(np.max(np.abs(s))),
568610 "setup_seconds": round(t1 - t0, 3),
611+ **({"strikes": strikes, "note": "experiment with a strike schedule, not a solution of the challenge"} if strikes else {}),
569612 }
570613 if fields is not None:
571614 import h5py