#!/usr/bin/env python3 """V1 validation ladder, steps 1-3 (areal-CMC gauge): flrw homogeneous FLRW to a chosen t/t_H: alpha = 1, A = 1, Kth = -H, particles static in x, P ∝ 1/a idcheck Yoo spherical CMC initial data vs the exact LTB functions (mass and energy per shell) ltb cold collapse for the Yoo profile compared shell-by-shell with the LTB cycloid in proper time ah apparent-horizon detection for mu = 0.5 (LTB: central singularity at 7.5 t_H) Usage: python3 scripts/v1_ladder.py {flrw,idcheck,ltb,ah} [--mu 0.3] [--tend 6] [--K 2000] [--n 40000] [--xout 8] """ import argparse, json, sys, time from pathlib import Path import numpy as np sys.path.insert(0, str(Path(__file__).resolve().parents[1])) from pbhgr.cosmo_ev import Background, CosmoGrid, CParticles, CosmoRun, CosmoRunConfig, solve_metric, metric_at from pbhgr.cosmo_id import yoo_spherical_cmc, flrw_profile, sample_shells, sample_shells_quiet from pbhgr.ltb import LTBGaussian ap = argparse.ArgumentParser() ap.add_argument("test"); ap.add_argument("--mu", type=float, default=0.3); ap.add_argument("--tend", type=float, default=6.0) ap.add_argument("--K", type=int, default=2000); ap.add_argument("--n", type=int, default=4); ap.add_argument("--xout", type=float, default=5.0) ap.add_argument("--cfl", type=float, default=0.4); ap.add_argument("--out", default="runs/v1"); ap.add_argument("--stretch", type=float, default=0.0); ap.add_argument("--sigma", type=float, default=0.0); ap.add_argument("--tinj", type=float, default=1.0); ap.add_argument("--xdisp", type=float, default=2.0); ap.add_argument("--diag", type=int, default=100); ap.add_argument("--stop_after_ah", type=float, default=0.0); ap.add_argument("--dR", type=float, default=0.0); ap.add_argument("--release", default=""); ap.add_argument("--rel_k", type=float, default=10.0); ap.add_argument("--rel_delta", type=float, default=0.05); ap.add_argument("--rel_h", type=float, default=0.6); ap.add_argument("--rel_xmax", type=float, default=2.0); ap.add_argument("--rel_bnl", type=float, default=0.0); ap.add_argument("--profile", default="yoo"); ap.add_argument("--n_inner", type=int, default=0); ap.add_argument("--x_inner", type=float, default=0.5) args = ap.parse_args() out = Path(args.out); out.mkdir(parents=True, exist_ok=True) bg = Background(H_i=5.0); rm = np.sqrt(6.0); t_H = bg.t_i * (np.sqrt(6) / 0.2) ** 3 from pbhgr.ltb import LTBProfile from pbhgr.profiles import profile_table _table = None if args.profile != "yoo": _table = profile_table(args.profile); _table["name"] = args.profile _Lref = LTBProfile(0.5 if args.test == "ah" else args.mu, _table); rm = _Lref.r_m; t_H = _Lref.t_H def make_L(mu_): return LTBGaussian(mu_) if _table is None else LTBProfile(mu_, _table) def make_prof(mu_, r_max): L_ = make_L(mu_) return yoo_spherical_cmc(mu_, r_max=r_max) if _table is None else yoo_spherical_cmc(mu_, r_max=r_max, zeta_fns=(L_.zeta, L_.dzeta, L_.ddzeta)) x_out = args.xout * rm grid = CosmoGrid(x_out=x_out, K=args.K, stretch=args.stretch) log = open(out / (f"{args.test}_mu{args.mu:.3f}_K{args.K}_n{args.n}" + (f"i{args.n_inner}" if args.n_inner else "") + f"_s{args.stretch:.0f}_sig{args.sigma:.0e}" + (f"_dR{args.dR:g}" if args.dR > 0 else "") + (f"_rel{args.release}_k{args.rel_k:g}_d{args.rel_delta:g}" + (f"_b{args.rel_bnl:g}" if args.rel_bnl > 0 else "") if args.release else "") + (f"_{args.profile}" if args.profile != "yoo" else "") + ".log"), "w") def P_(*a): print(*a); print(*a, file=log); log.flush() P_(f"t_i={bg.t_i:.5f} t_H={t_H:.3f} ({t_H/bg.t_i:.1f} t_i) x_out={x_out:.3f} (H_i x_out={5*x_out:.0f}) K={args.K} stretch={args.stretch} dx_min={grid.dx:.2e} ({grid.dx/rm:.1e} r_m) dx_max={grid.dxc.max():.2e} particles/cell={args.n}") def cb(run, d): if args.stop_after_ah > 0 and run.ah_first is not None and run.t > run.ah_first["t"] + args.stop_after_ah * t_H: run.cfg.t_end = run.t # classification done: stop shortly after the first apparent horizon ah = d["AH"] P_(f"step {d['step']:6d} t/t_H={d['t']/t_H:8.4f} a={d['a']:9.3f} s={d.get('stretch',0):.2f} dRmin={d.get('dR_min',0):.2f} max2m/R(coll)={d['max_2m_over_R_collapsing']:.4f} at x/rm={d['x_at_max']/rm:.3f} " f"min_alpha={d['min_alpha']:.4f} max_A={d['max_A']:.4f} m_out/flrw-1={d['m_out_over_flrw']-1:+.2e} bad={d['n_bad']} " f"ueq={d['unused_eq_residual']:.2e} min_x={d['min_x']:.2e} xbad/rm={d['x_bad_max']/rm:.4f}" + (f" AH: x/rm={ah['x_AH']/rm:.4f} M={ah['M_AH']:.4g}" if ah else "")) if args.test == "flrw": prof = flrw_profile(H_i=5.0, r_max=x_out * 1.01) x, P, Lsq, N = sample_shells_quiet(prof, grid.x_edge, args.n, n_inner=args.n_inner, x_inner=args.x_inner * rm) # quiet start, P = 0 pec = np.zeros(len(x), bool) part = CParticles(x, P, Lsq, N) met0 = solve_metric(grid, bg, part, bg.t_i) P_(f"t_i: max|alpha-1|={np.max(np.abs(met0.alpha_c-1)):.2e} max|A-1|={np.max(np.abs(met0.A_edge-1)):.2e} max|Kth+H|/H={np.max(np.abs(met0.Kth_edge+5))/5:.2e} m_out/flrw-1={met0.m_edge[-1]/(0.5*25*x_out**3)-1:+.2e}") cfg = CosmoRunConfig(cfl=args.cfl, t_end=args.tend * t_H, diag_every=200) run = CosmoRun(grid, bg, part, cfg, bg.t_i); x0 = part.x.copy(); P0 = part.P.copy() t0 = time.time(); met = run.run(cb) a = met.a P_(f"done: steps={run.step} wall={time.time()-t0:.1f}s max|alpha-1|={np.max(np.abs(met.alpha_c-1)):.2e} max|A-1|={np.max(np.abs(met.A_edge-1)):.2e} " f"max|Kth+H|/H={np.max(np.abs(met.Kth_edge+bg.H(run.t)))/bg.H(run.t):.2e} max|x-x0|/x_out={np.max(np.abs(part.x-x0))/x_out:.2e} " f"max|P| (rest)={np.max(np.abs(part.P)):.2e} m_out/flrw-1={met.m_edge[-1]/(0.5*bg.H(run.t)**2*met.R_edge[-1]**3)-1:+.2e}") elif args.test == "idcheck": prof = make_prof(args.mu, x_out * 1.05) L = make_L(args.mu) # enclosed rest mass as the shell label: LTB M_rest(r) = int dm/sqrt(1+2E); ID M_rest(R) = int dN/dR dR r = prof["r"]; R = prof["R"] Mrest_id = np.concatenate([[0], np.cumsum(0.5 * (prof["dN_dR"][1:] + prof["dN_dR"][:-1]) * np.diff(R))]) rL = np.linspace(1e-6, x_out * 1.05, 80001); mL = L.m(rL); EL = L.E(rL) Mrest_L = np.concatenate([[0], np.cumsum(np.diff(mL) / np.sqrt(1 + 2 * 0.5 * (EL[1:] + EL[:-1])))]) # shell energy from the ID state: dR/dtau = P/A - R Kth W with P = Gamma v, W = Gamma Gam = prof["Gamma"]; Pn = Gam * prof["v"] dRdtau = Pn / prof["A"] - R * prof["Kth"] * Gam E_id = 0.5 * dRdtau**2 - prof["m"] / np.maximum(R, 1e-300) tau0 = LTBGaussian.shell_tau(np.maximum(prof["m"], 1e-300), np.minimum(E_id, -1e-300), R, dRdtau) P_("Yoo CMC data vs exact LTB at equal enclosed rest mass (eps_i = 0.2):") P_(" r_iso R m_id/m_LTB-1 E_id/E_LTB-1 delta=E/rho_bar-1 v tau0/t_i-1") rho_bar = 3 * 25 / (8 * np.pi) for rr in (0.25, 0.5, 1.0, 1.5, 2.0, 2.449, 3.0, 4.0, 6.0, 8.0, 12.0): i = np.searchsorted(r, rr) Mi = Mrest_id[i] rl = np.interp(Mi, Mrest_L, rL) P_(f" {r[i]:6.3f} {R[i]:7.4f} {prof['m'][i]/L.m(rl)-1:+.3e} {E_id[i]/L.E(rl)-1:+.3e} {prof['E'][i]/rho_bar-1:+.4e} {prof['v'][i]:+.2e} {tau0[i]/bg.t_i-1:+.3e}") P_(f" total mass excess inside x_out relative to FLRW: {(prof['m'][np.searchsorted(r, x_out)]/(0.5*25*R[np.searchsorted(r, x_out)]**3)-1):+.3e}") P_(f" central delta = {prof['E'][0]/rho_bar-1:.5f} (linear growing mode: (2/5) eps^2 * 6/... = 2.4 mu eps_i^2/6*...; expected ~ (2/5)(k^2/(aH)^2)*mu = {0.4*0.04*args.mu:.5f})") elif args.test in ("ltb", "ah"): mu = 0.5 if args.test == "ah" else args.mu prof = make_prof(mu, x_out * 1.05) x, P, Lsq, N = sample_shells_quiet(prof, grid.x_edge, args.n, n_inner=args.n_inner, x_inner=args.x_inner * rm) part = CParticles(x, P, Lsq, N) L = make_L(mu) tC0 = L.tC(1e-6) if not np.isfinite(tC0) or tC0 <= 0: tC0 = L.tC(rm) # flat-core profiles: no central growing mode, use the r_m shell P_(f"mu={mu}: LTB t_C(0)={tC0/t_H:.3f} t_H, t_C(r_m)={L.tC(rm)/t_H:.2f} t_H; run to {args.tend} t_H") met0 = solve_metric(grid, bg, part, bg.t_i) R0 = met0.a * part.x # shell labels: conserved rest mass M_i = sum_{j 0 else 0.0, inject_sigma=args.sigma, inject_xmax=args.xdisp * rm, dR_target=args.dR, release_mode=args.release, release_k_ratio=args.rel_k, release_delta_ent=args.rel_delta, release_h=args.rel_h, release_t_ent=t_H, release_xmax=args.rel_xmax * rm, release_b_nl=args.rel_bnl) run = CosmoRun(grid, bg, part, cfg, bg.t_i) t0 = time.time(); met = run.run(cb) if args.release: P_(f" release mode {args.release}: k ratio {args.rel_k}, delta/zeta {args.rel_delta}, h {args.rel_h}; shells released {run.n_released}, first release at t/t_H={run.t_first_release/t_H:.2f}") if run.injected: P_(f" dispersion injected at t/t_H={run.injected['t']/t_H:.2f}: sigma={run.injected['sigma']:.1e} on {run.injected['n']} shells (x < {args.xdisp} r_m)") P_(f"done: steps={run.step} wall={time.time()-t0:.1f}s; first AH: {run.ah_first}") tr = ~np.isnan(run.tau_trap) if tr.any(): tt = tau_off + run.tau_trap tL = np.array([L.t_AH(r) for r in r_lab]) P_(f" per-shell trapping (proper time): {tr.sum()} shells trapped; labels r/r_m in [{r_lab[tr].min()/rm:.3f},{r_lab[tr].max()/rm:.3f}]; " f"earliest tau_trap={np.nanmin(tt[tr])/t_H:.2f} t_H (label r/r_m={r_lab[tr][np.argmin(tt[tr])]/rm:.3f}); median tau_trap/LTB t_AH = {np.median(tt[tr]/tL[tr]):.4f}, p5-p95 [{np.percentile(tt[tr]/tL[tr],5):.4f},{np.percentile(tt[tr]/tL[tr],95):.4f}]; " f"median coordinate lag t_detect - tau_trap = {np.median(run.t_trap[tr]-tt[tr])/t_H:.2f} t_H") for lo, hi in ((0,0.1),(0.1,0.3),(0.3,0.5),(0.5,0.8),(0.8,1.2)): sel = tr & (r_lab >= lo*rm) & (r_lab < hi*rm) if sel.any(): P_(f" labels [{lo},{hi}) r_m: n={sel.sum()}, tau_trap/t_AH^LTB median={np.median(tt[sel]/tL[sel]):.4f}, LTB t_AH range [{tL[sel].min()/t_H:.2f},{tL[sel].max()/t_H:.2f}] t_H") if run.ah_first is not None: sA = run.ah_snapshot; xA = run.ah_first["x_AH"] near = np.argsort(np.abs(sA["R"] / sA["a"] - xA))[:8] tauA = (tau_off + sA["tau"])[near] lab_AH = np.median(r_lab[near]); inside = r_lab < lab_AH; tauI = tau_off + sA["tau"] P_(f" shells inside the first AH (labels < {lab_AH/rm:.3f} r_m): n={inside.sum()}, proper time median={np.median(tauI[inside])/t_H:.2f} t_H, min={tauI[inside].min()/t_H:.2f}, max={tauI[inside].max()/t_H:.2f}; " f"fraction re-expanded beyond x_AH: {np.mean(sA['R'][inside]/sA['a'] > xA):.3f}; LTB caustic t_C(0)={L.tC(1e-6)/t_H:.2f}") P_(f" AH invariant check: coordinate t={run.ah_first['t']/t_H:.2f} t_H, M_AH={run.ah_first['M_AH']:.2f} (= {run.ah_first['M_AH']/(0.75*t_H):.3f} M_H(t_k)); shells at the AH: " f"label r/r_m={np.median(r_lab[near])/rm:.3f}, LTB m={np.median(m0[near]):.2f}, shell proper time={np.median(tauA)/t_H:.2f} t_H vs LTB t_AH(shell)={np.median([L.t_AH(r) for r in r_lab[near]])/t_H:.2f} t_H, LTB t_C(shell)={np.median([L.tC(r) for r in r_lab[near]])/t_H:.2f}") # LTB comparison at snapshots: predicted R(tau) from the shell invariants (m0, E0) fixed at t_i inner = (x < 1.5 * rm) & (E0 < -1e-10) & (m0 > 0) & (r_lab > 0.02) rows = [] for s in run.snapshots: tau = (tau_off + s["tau"])[inner] Rpred = LTBGaussian.shell_R_of_tau(m0[inner], E0[inner], tau) err = np.abs(s["R"][inner] / Rpred - 1) mdrift = np.abs(s["m"][inner] / m0[inner] - 1) inner_ = np.ones(inner.sum(), bool) rows.append((s["t"] / t_H, s["t"] / tC0, np.median(err), np.percentile(err, 95), np.median(mdrift), np.percentile(mdrift, 95))) P_(f" snapshot t/t_H={rows[-1][0]:7.3f} (t/t_C(0)={rows[-1][1]:.2f}): |R/R_LTB-1| median={rows[-1][2]:.2e} p95={rows[-1][3]:.2e}; enclosed-mass drift median={rows[-1][4]:.2e} p95={rows[-1][5]:.2e} (shells x<1.5 r_m)") s = run.snapshots[-1]; tau = (tau_off + s["tau"])[inner] Rpred = LTBGaussian.shell_R_of_tau(m0[inner], E0[inner], tau); err = s["R"][inner] / Rpred - 1; xi = x[inner] P_(" last snapshot, error vs radius (bins in x/r_m): median |R/R_LTB-1|, p95, median m-drift") for lo, hi in ((0, 0.05), (0.05, 0.2), (0.2, 0.5), (0.5, 1.0), (1.0, 1.5)): sel = (xi >= lo * rm) & (xi < hi * rm) if sel.any(): P_(f" [{lo:.2f},{hi:.2f}) n={sel.sum():5d}: {np.median(np.abs(err[sel])):.2e} {np.percentile(np.abs(err[sel]),95):.2e} {np.median(np.abs(s['m'][inner][sel]/m0[inner][sel]-1)):.2e} E0 in [{E0[inner][sel].min():.3f},{E0[inner][sel].max():.3f}]") if hasattr(run, "ueq_profile"): up = np.abs(run.ueq_profile); xe = grid.x_edge P_(" unused-equation residual profile (|res|/scale) by x/r_m: " + ", ".join(f"[{lo},{hi}):{np.max(up[(xe>=lo*rm)&(xe 0 else "") + (f"_rel{args.release}_k{args.rel_k:g}_d{args.rel_delta:g}" + (f"_b{args.rel_bnl:g}" if args.rel_bnl > 0 else "") if args.release else "") + (f"_{args.profile}" if args.profile != "yoo" else "") json.dump({"diag": run.diag, "ltb_rows": rows, "ah_first": run.ah_first}, open(out / f"{tag}.json", "w"), default=float) np.savez(out / f"{tag}_shells.npz", x0=x, N=N, r_lab=r_lab, m0=m0, E0=E0, tau_off=tau_off, tau_trap=run.tau_trap, t_trap=run.t_trap, tau_end=part.tau, x_end=part.x, P_end=part.P, t_H=t_H, mu=mu)