70cc577Struck FPUT chain: model, experiments, figures, web appJeremy Magland 1"""Experiments with the struck FPUT-beta chain. Each writes results/<name>.npz; plots.py draws them.
3 python experiments.py scan # Lyapunov exponent and equipartition vs strike energy (undamped)
4 python experiments.py cascade # mode energies over time at three strike energies (undamped)
5 python experiments.py damped # damped strikes: ensembles of 16 nearby strikes, frozen outcome
6 python experiments.py all
8N = 32, beta = 1, strike centered on mass 8 with width 2 masses. mu = 0 (harmonic at low k) or mu = 1
9(on-site springs: an inharmonic spectrum, closer to the plate).
10"""
12from __future__ import annotations
14import os
15import sys
16from multiprocessing import Pool
18os.environ.setdefault("OMP_NUM_THREADS", "1")
20import numpy as np # noqa: E402
22from model import Chain, run, spectral_entropy, strike # noqa: E402
24HERE = os.path.dirname(os.path.abspath(__file__))
25RES = os.path.join(HERE, "results")
26N, CENTER, WIDTH = 32, 8, 2.0
29def dt_for(E, mu):
30 return min(0.1, 0.25 / np.sqrt(4 + mu**2 + 12 * E))
33def _scan_one(args):
34 mu, E = args
35 ch = Chain(N=N, mu=mu)
36 M = 6
37 u0, v0 = strike(ch, E, CENTER, WIDTH, M=M, scale=1 + 1e-6 * np.arange(M))
38 dt = dt_for(E, mu)
39 r = run(ch, u0, v0, 4000.0, dt, lyap=True, record_every=int(round(20 / dt)))
40 neff = spectral_entropy(r["Ek"]) # (T, M)
41 return mu, E, r["lyap"], r["t"], neff
44def scan(pool):
45 energies = np.geomspace(0.05, 50, 16)
46 out = {}
47 for mu, E, lam, t, neff in pool.imap_unordered(_scan_one, [(mu, E) for mu in (0.0, 1.0) for E in energies]):
48 i = int(np.argmin(abs(energies - E)))
49 out[f"mu{mu:g}_E{i}_lyap"], out[f"mu{mu:g}_E{i}_neff"], out[f"mu{mu:g}_t"] = lam, neff, t
50 print(f" scan mu {mu:g} E {E:.3g}: lambda {np.median(lam):.4f}, n_eff/N final {np.median(neff[-1]):.2f}", flush=True)
51 np.savez(os.path.join(RES, "scan.npz"), energies=energies, **out)
54def _cascade_one(args):
55 mu, E = args
56 ch = Chain(N=N, mu=mu)
57 u0, v0 = strike(ch, E, CENTER, WIDTH)
58 dt = dt_for(E, mu)
59 r = run(ch, u0, v0, 3000.0, dt, record_every=int(round(2 / dt)))
60 return mu, E, r["t"], r["Ek"][:, 0, :], ch.omega
63def cascade(pool):
64 out = {}
65 jobs = [(mu, E) for mu in (0.0, 1.0) for E in (0.5, 2.0, 10.0)]
66 for mu, E, t, Ek, om in pool.imap_unordered(_cascade_one, jobs):
67 k = f"mu{mu:g}_E{E:g}"
68 out[k + "_t"], out[k + "_Ek"], out[k + "_omega"] = t, Ek, om
69 print(f" cascade mu {mu:g} E {E:g}", flush=True)
70 np.savez(os.path.join(RES, "cascade.npz"), **out)
73def damping_for(ch, sigma_low, sigma_high):
74 """c0, c1 so that the lowest mode decays (amplitude) at sigma_low and the highest at sigma_high."""
75 om = ch.omega
76 a1, aN = om[0] ** 2 - ch.mu**2, om[-1] ** 2 - ch.mu**2
77 c1 = 2 * (sigma_high - sigma_low) / (aN - a1)
78 c0 = 2 * sigma_low - c1 * a1
79 return c0, c1
82def _damped_one(args):
83 mu, E0 = args
84 base = Chain(N=N, mu=mu)
85 c0, c1 = damping_for(base, 4e-4, 3e-3)
86 ch = Chain(N=N, mu=mu, c0=c0, c1=c1)
87 M = 16
88 u0, v0 = strike(ch, E0, CENTER, WIDTH, M=M, scale=1 + 1e-10 * np.arange(M))
89 dt = dt_for(E0, mu)
90 r = run(ch, u0, v0, 8000.0, dt, lyap=True, record_every=int(round(5 / dt)))
91 return mu, E0, r["t"], r["Ek"], r["E"], r["logsum"], c0, c1
94def damped(pool):
95 out = {}
96 jobs = [(mu, E0) for mu in (0.0, 1.0) for E0 in (0.5, 5.0, 30.0)]
97 for mu, E0, t, Ek, E, logsum, c0, c1 in pool.imap_unordered(_damped_one, jobs):
98 k = f"mu{mu:g}_E{E0:g}"
99 out[k + "_t"], out[k + "_Ek"], out[k + "_E"], out[k + "_logsum"] = t, Ek, E, logsum
100 out[k + "_c"] = np.array([c0, c1])
101 print(f" damped mu {mu:g} E0 {E0:g}", flush=True)
102 np.savez(os.path.join(RES, "damped.npz"), **out)
105if __name__ == "__main__":
106 os.makedirs(RES, exist_ok=True)
107 which = sys.argv[1:] or ["all"]
108 if "all" in which:
109 which = ["scan", "cascade", "damped"]
110 with Pool(int(os.environ.get("WORKERS", "16"))) as pool:
111 for w in which:
112 print(w, flush=True)
113 dict(scan=scan, cascade=cascade, damped=damped)[w](pool)