diff --git a/analysis/BENCHMARK_NOTES.md b/analysis/BENCHMARK_NOTES.md index ce98be0..ce45a38 100644 --- a/analysis/BENCHMARK_NOTES.md +++ b/analysis/BENCHMARK_NOTES.md @@ -29,13 +29,13 @@ single source of truth; the sweep script and the per-flap configs must agree. ## Hydrodynamics - **Physical plant coordinate:** the time-domain solver measures `P_capture` in the **hinge DOF**. -- **Current post-processing basis:** the committed benchmark still tabulates - `F_exc_Nm` / `B55_Nmsrad` from the CG-referenced H5 files - (`hydroData/vgoswec_{0,10,20,45,90}.h5`). -- **Implication:** `P_opt = |F_exc|^2 / (8 * B55)` is invariant under the CG↔hinge - referral, so the efficiency denominator is only mildly affected, but the - **masking** rule is basis-sensitive. A hinge-basis re-tabulation remains an open - post-processing task. +- **Post-processing basis:** the committed capture-efficiency CSVs now tabulate + `F_exc_Nm` / `B55_Nmsrad` from the **hinge-referenced** H5 files + (`hydroData/hinged_vgoswec_{0,10,20,45,90}.h5`), matching the hinge DOF used by + the controller and `P_capture`. +- **Implication:** `P_opt = |F_exc|^2 / (8 * B55)` remains basis-invariant in exact + referral, and the `masked` rule is now applied on the same hinge basis as the + time-domain capture curves. - **Hinge-referenced impedance files now matter elsewhere:** controller gain computation (`opt_passive`, CC) now explicitly uses `hydro.impedance_h5_file = hydroData/hinged_vgoswec_*.h5`. See @@ -43,9 +43,10 @@ single source of truth; the sweep script and the per-flap configs must agree. - **P_opt (Budal / Falnes bound):** `P_opt = |F_exc|^2 / (8 * B55)`, computed from body1 pitch hydro (`radiation_damping/components/5_5` and `excitation/mag[dof=5,dir=0]`), at H = 0.05 m. -- **De-normalization (WEC-Sim / BEMIO convention, rho and g read from each H5, - rho = 1000, g = 9.80665):** - - `B55 = B55_norm * rho * omega` [N*m/(rad/s)] (peak ~3, matches BEMRosetta) +- **De-normalization (WEC-Sim / BEMIO convention, `rho` and `g` read from each H5; + hinged files currently store `rho = 1025`, whereas the earlier CG-basis notes + assumed `rho = 1000`):** + - `B55 = max(0, B55_norm * rho * omega)` [N*m/(rad/s)] (peak ~3, matches BEMRosetta) - `F_exc = F_exc_norm * rho * g * a` [N*m] (per-amplitude ~174 N*m/m, matches BEMRosetta) - **Mask rule:** points with `B55 <= 1e-04 N*m*s/rad` are treated as reactive-limited and omitted. On the 4-12 band all five flaps clear this @@ -69,15 +70,11 @@ single source of truth; the sweep script and the per-flap configs must agree. omega = 7.0 and 7.5). At the 0.5 rad/s grid spacing this edge is under-resolved, so the exact peak height (~34 %) is grid-sensitive; the shape and location are robust. Refining the grid near omega ~ 7 would pin the peak value. -- **CG- vs hinge-referenced:** P_opt/eta here use CG-referenced free-flap pitch - hydro, not the hinge-referenced flap. `F_exc` is the raw un-hinge-referred pitch - moment and `alpha` is signed positive to paper over the resulting phase/sign - mismatch. +- **Hinge-referenced tabulation:** P_opt/eta and the `masked` rule are now tabulated + from hinge-referenced flap hydro, so the benchmark notes no longer rely on a + CG-basis workaround or sign-convention patch-up in post-processing. ## Follow-up (next milestone) -- Re-tabulate `F_exc`, `B55`, and the `masked` column from the - **hinge-referenced** coefficients (`hinged_vgoswec_*.h5`) so benchmark masking is - consistent with the hinge DOF used by the plant. - Implement true complex-conjugate control (theoretical eta_max reference) with hinge-referenced K_r / B_r, then causal approximations (Korde) and constrained MPC (Ringwood). Resonance markers return once resonance is defined consistently diff --git a/analysis/README.md b/analysis/README.md index 7c80ed2..406a000 100644 --- a/analysis/README.md +++ b/analysis/README.md @@ -57,10 +57,10 @@ | `omega_rads` | Angular frequency [rad/s], `2π/T` | | `P_capture_W` | Steady-state mean absorbed power from tuned `exc_ff_pid` [W] | | `P_opt_W` | Theoretical optimum power [W], blank when masked | -| `B55_Nmsrad` | De-normalized pitch radiation damping `B55` [N·m·s/rad] | -| `F_exc_Nm` | De-normalized pitch excitation moment magnitude `|F_exc|` for `A=0.025 m` [N·m] | +| `B55_Nmsrad` | Hinge-basis de-normalized pitch radiation damping `B55` [N·m·s/rad] from `hydroData/hinged_vgoswec_*.h5` | +| `F_exc_Nm` | Hinge-basis de-normalized pitch excitation moment magnitude `|F_exc|` for `A=0.025 m` [N·m] from `hydroData/hinged_vgoswec_*.h5` | | `eta` | Capture efficiency `η = P_capture / P_opt`; only blank when `masked=true` (`B55 <= 1e-4`) | -| `masked` | `true` where `B55 <= 1e-4` (reactive-limited / undefined `P_opt`, including non-positive `B55`) | +| `masked` | `true` where hinge-basis `B55 <= 1e-4` (reactive-limited / undefined `P_opt`, including non-positive `B55`) | ## Capture-efficiency method (tuned `exc_ff_pid`) @@ -75,7 +75,9 @@ - Body: `body1` (flap), ignore `body2`. - Pitch term: `body1/hydro_coeffs/radiation_damping/components/5_5`. - Excitation: `body1/hydro_coeffs/excitation/mag` at DOF5 (index 4), direction 0. - - De-normalization: `B55 = B55_norm * rho * omega`, `|F_exc| = mag * rho * g * A`. + - Basis: hinge-referenced `hydroData/hinged_vgoswec_*.h5`, matching the hinge DOF used by `P_capture_W`. + - De-normalization: `B55 = max(0, B55_norm * rho * omega)`, `|F_exc| = mag * rho * g * A`. + - `rho` and `g` are read from each H5 file; the hinged files store `rho = 1025`. - Wave amplitude fixed to `A = 0.025 m` (`H = 0.05 m`) for both sim and `P_opt`. - Masking/flagging (essential): - `B55 <= 1e-4` => `P_opt` undefined (reactive-limited notch), so `η` is not reported/plotted (`masked=true`). @@ -112,7 +114,10 @@ python3 scripts/plot_kpkd_surface.py python3 scripts/capture_efficiency_sweep.py # 4. Re-plot only, from committed capture-efficiency CSVs +python3 scripts/retabulate_hydro_columns.py --dry-run +python3 scripts/retabulate_hydro_columns.py python3 scripts/capture_efficiency_sweep.py --plot-only python3 scripts/cc_capture_efficiency_sweep.py --plot-only python3 scripts/cc_vs_ffpid_comparison.py --plot-only +python3 scripts/plot_opt_passive_stage.py ``` diff --git a/docs/EOD_SUMMARY_2026-09-19.md b/docs/EOD_SUMMARY_2026-09-19.md index 57dd90c..b3a8b20 100644 --- a/docs/EOD_SUMMARY_2026-09-19.md +++ b/docs/EOD_SUMMARY_2026-09-19.md @@ -131,10 +131,16 @@ The owner should be able to `git pull` and immediately re-run the `opt_passive` 1. **Highest priority:** ff+PID guard-fire fraction was never measured. Instrument `ExcitationVelocityController::ComputeForce` to count `tau*vel > 0` events and run VGM-0 at `T = 4.50 s`. If the guard fires near 100% of the time, the long-period band is the passive-safety floor rather than feedforward control, and the claim *"ff+PID carries the long-period tail"* must be rewritten. 2. The `2.44488972 W` headline is still unconfirmed as settled (`430 s` vs `760 s` check still needed). -3. Sweep scripts still tabulate `F_exc_Nm` / `B55_Nmsrad` from **CG** H5 files while `P_capture` is measured in the hinge DOF. `P_opt = F^2/(8B)` is basis-invariant so `eta` is only mildly affected, but `masked` is not. Hinge-basis re-tabulation is still needed as post-processing. +3. **Resolved in follow-up PR:** sweep post-processing now reads hinge-referenced + `hydroData/hinged_vgoswec_*.h5`, clamps de-normalized `B55 >= 0` in Python to + match `src/impedance.cpp`, and can re-tabulate existing `analysis/{passive,opt_passive}` + CSVs in place without re-running simulations (`scripts/retabulate_hydro_columns.py`). 4. `analysis/FINDINGS_3REGIME.md` hand-written tables are still stale. 5. There remains a uniform `~1.4%` residual between corrected analytic `omega_n` and free-decay (`single-DOF analytic model` vs `coupled plant`). 6. `[impedance] INFO: legacy A55-match rho=...` still prints more often than intended and should become truly once-per-process. +7. The isolated single-point `P_capture` dips near resonance (notably VGM-0 `T=6.00 s`, + VGM-45 `T=3.50 s`, VGM-90 `T=3.00 s`) remain an open question; no controller or + hydro post-processing change in this PR attempts to explain or alter them. --- diff --git a/docs/img/opt_passive_stage_all_flaps.png b/docs/img/opt_passive_stage_all_flaps.png new file mode 100644 index 0000000..2af4c97 Binary files /dev/null and b/docs/img/opt_passive_stage_all_flaps.png differ diff --git a/docs/img/opt_passive_stage_gain.png b/docs/img/opt_passive_stage_gain.png new file mode 100644 index 0000000..1c51de2 Binary files /dev/null and b/docs/img/opt_passive_stage_gain.png differ diff --git a/docs/img/opt_passive_stage_per_flap.png b/docs/img/opt_passive_stage_per_flap.png new file mode 100644 index 0000000..1f40cff Binary files /dev/null and b/docs/img/opt_passive_stage_per_flap.png differ diff --git a/scripts/capture_efficiency_sweep.py b/scripts/capture_efficiency_sweep.py index 9747085..234613b 100755 --- a/scripts/capture_efficiency_sweep.py +++ b/scripts/capture_efficiency_sweep.py @@ -65,27 +65,27 @@ 0: { "label": "VGM-0", "config": "config/vgoswec_0_exc_ff_pid.yaml", - "h5": "hydroData/vgoswec_0.h5", + "h5": "hydroData/hinged_vgoswec_0.h5", }, 10: { "label": "VGM-10", "config": "config/vgoswec_10_exc_ff_pid.yaml", - "h5": "hydroData/vgoswec_10.h5", + "h5": "hydroData/hinged_vgoswec_10.h5", }, 20: { "label": "VGM-20", "config": "config/vgoswec_20_exc_ff_pid.yaml", - "h5": "hydroData/vgoswec_20.h5", + "h5": "hydroData/hinged_vgoswec_20.h5", }, 45: { "label": "VGM-45", "config": "config/vgoswec_45_exc_ff_pid.yaml", - "h5": "hydroData/vgoswec_45.h5", + "h5": "hydroData/hinged_vgoswec_45.h5", }, 90: { "label": "VGM-90", "config": "config/vgoswec_90_exc_ff_pid.yaml", - "h5": "hydroData/vgoswec_90.h5", + "h5": "hydroData/hinged_vgoswec_90.h5", }, } @@ -208,7 +208,7 @@ def popt_curve_from_h5(h5_path: Path, periods_s: np.ndarray) -> tuple[np.ndarray b55_norm = b55_norm[order] fexc_norm = fexc_norm[order] - b55 = b55_norm * rho * w + b55 = np.maximum(0.0, b55_norm * rho * w) fexc = fexc_norm * rho * g * WAVE_AMPLITUDE_M omega_targets = (2.0 * math.pi) / periods_s diff --git a/scripts/cc_capture_efficiency_sweep.py b/scripts/cc_capture_efficiency_sweep.py index 35790d7..77c8360 100644 --- a/scripts/cc_capture_efficiency_sweep.py +++ b/scripts/cc_capture_efficiency_sweep.py @@ -55,11 +55,11 @@ REACTIVE_NOTE = f"|P_capture| / P_converted < {REACTIVE_RESIDUAL_TOL:.0e}" FLAPS = { - 0: {"label": "VGM-0", "config": "config/vgoswec_0_cc.yaml", "h5": "hydroData/vgoswec_0.h5"}, - 10: {"label": "VGM-10", "config": "config/vgoswec_10_cc.yaml", "h5": "hydroData/vgoswec_10.h5"}, - 20: {"label": "VGM-20", "config": "config/vgoswec_20_cc.yaml", "h5": "hydroData/vgoswec_20.h5"}, - 45: {"label": "VGM-45", "config": "config/vgoswec_45_cc.yaml", "h5": "hydroData/vgoswec_45.h5"}, - 90: {"label": "VGM-90", "config": "config/vgoswec_90_cc.yaml", "h5": "hydroData/vgoswec_90.h5"}, + 0: {"label": "VGM-0", "config": "config/vgoswec_0_cc.yaml", "h5": "hydroData/hinged_vgoswec_0.h5"}, + 10: {"label": "VGM-10", "config": "config/vgoswec_10_cc.yaml", "h5": "hydroData/hinged_vgoswec_10.h5"}, + 20: {"label": "VGM-20", "config": "config/vgoswec_20_cc.yaml", "h5": "hydroData/hinged_vgoswec_20.h5"}, + 45: {"label": "VGM-45", "config": "config/vgoswec_45_cc.yaml", "h5": "hydroData/hinged_vgoswec_45.h5"}, + 90: {"label": "VGM-90", "config": "config/vgoswec_90_cc.yaml", "h5": "hydroData/hinged_vgoswec_90.h5"}, } JOURNAL_STYLE = { @@ -199,7 +199,7 @@ def popt_curve_from_h5(h5_path: Path, periods_s: np.ndarray) -> tuple[np.ndarray omega_rads = omega_rads[order] b55_norm = b55_norm[order] fexc_norm = fexc_norm[order] - b55 = b55_norm * rho * omega_rads + b55 = np.maximum(0.0, b55_norm * rho * omega_rads) fexc = fexc_norm * rho * g * WAVE_AMPLITUDE_M omega_targets = (2.0 * math.pi) / periods_s diff --git a/scripts/passive_vs_optpassive_sweep.py b/scripts/passive_vs_optpassive_sweep.py index dca3306..095a35a 100644 --- a/scripts/passive_vs_optpassive_sweep.py +++ b/scripts/passive_vs_optpassive_sweep.py @@ -83,35 +83,35 @@ "label": "VGM-0", "passive_config": "config/vgoswec_0_passive.yaml", "opt_passive_config": "config/vgoswec_0_opt_passive.yaml", - "h5": "hydroData/vgoswec_0.h5", + "h5": "hydroData/hinged_vgoswec_0.h5", "omega0": 1.07, }, 10: { "label": "VGM-10", "passive_config": "config/vgoswec_10_passive.yaml", "opt_passive_config": "config/vgoswec_10_opt_passive.yaml", - "h5": "hydroData/vgoswec_10.h5", + "h5": "hydroData/hinged_vgoswec_10.h5", "omega0": 1.468, }, 20: { "label": "VGM-20", "passive_config": "config/vgoswec_20_passive.yaml", "opt_passive_config": "config/vgoswec_20_opt_passive.yaml", - "h5": "hydroData/vgoswec_20.h5", + "h5": "hydroData/hinged_vgoswec_20.h5", "omega0": 1.568, }, 45: { "label": "VGM-45", "passive_config": "config/vgoswec_45_passive.yaml", "opt_passive_config": "config/vgoswec_45_opt_passive.yaml", - "h5": "hydroData/vgoswec_45.h5", + "h5": "hydroData/hinged_vgoswec_45.h5", "omega0": 1.84, }, 90: { "label": "VGM-90", "passive_config": "config/vgoswec_90_passive.yaml", "opt_passive_config": "config/vgoswec_90_opt_passive.yaml", - "h5": "hydroData/vgoswec_90.h5", + "h5": "hydroData/hinged_vgoswec_90.h5", "omega0": 2.094, }, } @@ -264,7 +264,7 @@ def popt_curve_from_h5( b55_norm = b55_norm[order] fexc_norm = fexc_norm[order] - b55 = b55_norm * rho * w + b55 = np.maximum(0.0, b55_norm * rho * w) fexc = fexc_norm * rho * g * WAVE_AMPLITUDE_M omega_targets = (2.0 * math.pi) / periods_s diff --git a/scripts/plot_opt_passive_stage.py b/scripts/plot_opt_passive_stage.py new file mode 100644 index 0000000..f67f8c1 --- /dev/null +++ b/scripts/plot_opt_passive_stage.py @@ -0,0 +1,189 @@ +#!/usr/bin/env python3 +"""Plot passive vs opt_passive stage figures from committed CSVs.""" + +from __future__ import annotations + +import argparse +import csv +import math +from pathlib import Path + +import matplotlib +matplotlib.use("Agg") +import matplotlib.pyplot as plt +import numpy as np + +from passive_vs_optpassive_sweep import FLAPS, JOURNAL_STYLE, load_efficiency_csv + + +def _load_resonance_periods(csv_path: Path) -> dict[int, float]: + with csv_path.open(newline="") as fh: + reader = csv.DictReader(fh) + if not reader.fieldnames or "angle_deg" not in reader.fieldnames or "cpp_zerocross_wn_rads" not in reader.fieldnames: + raise RuntimeError( + f"Expected columns 'angle_deg' and 'cpp_zerocross_wn_rads' in {csv_path}" + ) + periods: dict[int, float] = {} + for row in reader: + wn = float(row["cpp_zerocross_wn_rads"]) + periods[int(row["angle_deg"])] = (2.0 * math.pi) / wn + if not periods: + raise RuntimeError(f"No resonance rows found in {csv_path}") + return periods + + +def _csv_map(repo: Path, subdir: str) -> dict[int, Path]: + return { + angle: repo / "analysis" / subdir / f"capture_efficiency_VGM{angle}.csv" + for angle in FLAPS + } + + +def _style_axes(ax) -> None: + ax.grid(True, which="major", linestyle="--", alpha=0.5) + ax.grid(True, which="minor", linestyle=":", alpha=0.35) + ax.minorticks_on() + ax.set_axisbelow(True) + + +def _paired_capture_arrays( + passive_csv: Path, + opt_csv: Path, +) -> tuple[np.ndarray, np.ndarray, np.ndarray]: + passive_rows = load_efficiency_csv(passive_csv) + opt_rows = load_efficiency_csv(opt_csv) + t_passive = np.array([row["T_s"] for row in passive_rows], dtype=float) + t_opt = np.array([row["T_s"] for row in opt_rows], dtype=float) + if t_passive.shape != t_opt.shape or not np.array_equal(t_passive, t_opt): + raise RuntimeError( + f"Mismatched T_s grids between {passive_csv.name} and {opt_csv.name}" + ) + p_passive = np.array([row["P_capture_W"] for row in passive_rows], dtype=float) + p_opt = np.array([row["P_capture_W"] for row in opt_rows], dtype=float) + return t_passive, p_passive, p_opt + + +def plot_per_flap( + passive_map: dict[int, Path], + opt_map: dict[int, Path], + t_res: dict[int, float], + out_png: Path, +) -> None: + fig, axes = plt.subplots(2, 3, figsize=(11.0, 6.5), sharex=True, sharey=True) + axes_flat = axes.flatten() + + for idx, angle in enumerate(sorted(FLAPS)): + ax = axes_flat[idx] + t_values, p_passive, p_opt = _paired_capture_arrays(passive_map[angle], opt_map[angle]) + + ax.plot(t_values, p_passive, linestyle="--", linewidth=1.8, color="tab:blue", label="passive") + ax.plot(t_values, p_opt, linestyle="-", linewidth=1.8, color="tab:orange", label="opt_passive") + ax.axvline(t_res[angle], linestyle=":", linewidth=1.4, color="0.25", label="$T_{res}$") + ax.set_title(FLAPS[angle]["label"]) + _style_axes(ax) + + axes_flat[0].legend(loc="best", fontsize=8) + for ax in axes[1, :]: + ax.set_xlabel("Wave period $T$ [s]") + for ax in axes[:, 0]: + ax.set_ylabel("$P_{capture}$ [W]") + axes_flat[-1].axis("off") + + fig.suptitle("Passive vs opt_passive capture power by flap") + fig.tight_layout() + out_png.parent.mkdir(parents=True, exist_ok=True) + fig.savefig(out_png) + plt.close(fig) + print(f"[ok] wrote {out_png}") + + +def plot_all_flaps(opt_map: dict[int, Path], t_res: dict[int, float], out_png: Path) -> None: + fig, ax = plt.subplots(figsize=(8.8, 4.8)) + colors = plt.cm.viridis(np.linspace(0.15, 0.9, len(FLAPS))) + + for color, angle in zip(colors, sorted(FLAPS)): + rows = load_efficiency_csv(opt_map[angle]) + t = np.array([row["T_s"] for row in rows], dtype=float) + p = np.array([row["P_capture_W"] for row in rows], dtype=float) + ax.plot(t, p, linewidth=1.9, color=color, label=FLAPS[angle]["label"]) + ax.axvline(t_res[angle], linestyle=":", linewidth=1.0, color=color, alpha=0.9) + + ax.set_xlabel("Wave period $T$ [s]") + ax.set_ylabel("$P_{capture}$ [W]") + ax.set_title("opt_passive capture power across all flaps") + _style_axes(ax) + ax.legend(loc="best", fontsize=8, ncol=2) + + fig.tight_layout() + out_png.parent.mkdir(parents=True, exist_ok=True) + fig.savefig(out_png) + plt.close(fig) + print(f"[ok] wrote {out_png}") + + +def plot_stage_gain(passive_map: dict[int, Path], opt_map: dict[int, Path], out_png: Path) -> None: + fig, ax = plt.subplots(figsize=(8.8, 4.8)) + colors = plt.cm.viridis(np.linspace(0.15, 0.9, len(FLAPS))) + + for color, angle in zip(colors, sorted(FLAPS)): + t, p_passive, p_opt = _paired_capture_arrays(passive_map[angle], opt_map[angle]) + ratio = np.divide( + p_opt, + p_passive, + out=np.full_like(p_opt, np.nan, dtype=float), + where=np.isfinite(p_opt) & np.isfinite(p_passive) & (p_passive > 0.0), + ) + valid = np.isfinite(ratio) & (ratio > 0.0) + ax.plot(t[valid], ratio[valid], linewidth=1.8, color=color, label=FLAPS[angle]["label"]) + + ax.axhline(1.0, linestyle="--", linewidth=1.2, color="0.2") + ax.set_xlabel("Wave period $T$ [s]") + ax.set_ylabel(r"$P_{\mathrm{opt\,passive}} / P_{\mathrm{passive}}$ [-]") + ax.set_title("opt_passive stage gain over passive") + ax.set_yscale("log") + _style_axes(ax) + ax.legend(loc="best", fontsize=8, ncol=2) + + fig.tight_layout() + out_png.parent.mkdir(parents=True, exist_ok=True) + fig.savefig(out_png) + plt.close(fig) + print(f"[ok] wrote {out_png}") + + +def parse_args() -> argparse.Namespace: + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument( + "--repo", + default=str(Path(__file__).resolve().parents[1]), + help="Repository root (default: parent of scripts/)", + ) + return parser.parse_args() + + +def main() -> int: + args = parse_args() + plt.rcParams.update(JOURNAL_STYLE) + repo = Path(args.repo).resolve() + + passive_map = _csv_map(repo, "passive") + opt_map = _csv_map(repo, "opt_passive") + missing = [path for path in [*passive_map.values(), *opt_map.values()] if not path.exists()] + if missing: + print(f"ERROR: Missing capture-efficiency CSVs, first missing file: {missing[0]}") + return 2 + + t_res = _load_resonance_periods(repo / "docs" / "freedecay_validation.csv") + missing_res = [angle for angle in FLAPS if angle not in t_res] + if missing_res: + print(f"ERROR: Missing resonance rows for flap angles: {missing_res}") + return 2 + out_dir = repo / "docs" / "img" + plot_per_flap(passive_map, opt_map, t_res, out_dir / "opt_passive_stage_per_flap.png") + plot_all_flaps(opt_map, t_res, out_dir / "opt_passive_stage_all_flaps.png") + plot_stage_gain(passive_map, opt_map, out_dir / "opt_passive_stage_gain.png") + return 0 + + +if __name__ == "__main__": + raise SystemExit(main()) diff --git a/scripts/retabulate_hydro_columns.py b/scripts/retabulate_hydro_columns.py new file mode 100644 index 0000000..2807898 --- /dev/null +++ b/scripts/retabulate_hydro_columns.py @@ -0,0 +1,141 @@ +#!/usr/bin/env python3 +"""Re-tabulate hydro-derived capture-efficiency columns without re-running sims.""" + +from __future__ import annotations + +import argparse +import csv +from collections.abc import Iterator +from pathlib import Path + +import numpy as np + +from passive_vs_optpassive_sweep import FLAPS, popt_curve_from_h5 + +TARGET_COLUMNS = ("B55_Nmsrad", "F_exc_Nm", "P_opt_W", "masked", "eta") +REQUIRED_COLUMNS = ("T_s", "P_capture_W", *TARGET_COLUMNS) +CSV_GROUPS = (("passive", "analysis/passive"), ("opt_passive", "analysis/opt_passive")) + + +def _format_float(value: float) -> str: + return f"{value:.8e}" + + +def _load_csv_exact(csv_path: Path) -> tuple[list[str], list[dict[str, str]]]: + with csv_path.open(newline="") as fh: + reader = csv.DictReader(fh) + rows = list(reader) + if not reader.fieldnames: + raise RuntimeError(f"CSV is missing a header row: {csv_path}") + if not rows: + raise RuntimeError(f"CSV has no rows: {csv_path}") + return list(reader.fieldnames), rows + + +def _retabulate_rows(rows: list[dict[str, str]], h5_path: Path) -> list[dict[str, str]]: + periods_s = np.array([float(row["T_s"]) for row in rows], dtype=float) + _, b55, fexc, p_opt, masked = popt_curve_from_h5(h5_path, periods_s) + + updated_rows: list[dict[str, str]] = [] + for idx, row in enumerate(rows): + updated = dict(row) + updated["B55_Nmsrad"] = _format_float(float(b55[idx])) + updated["F_exc_Nm"] = _format_float(float(fexc[idx])) + updated["masked"] = "true" if bool(masked[idx]) else "false" + if bool(masked[idx]): + updated["P_opt_W"] = "" + updated["eta"] = "" + else: + updated["P_opt_W"] = _format_float(float(p_opt[idx])) + p_capture_text = row.get("P_capture_W", "").strip() + if p_capture_text: + p_capture = float(p_capture_text) + if np.isfinite(p_capture) and np.isfinite(p_opt[idx]) and float(p_opt[idx]) > 0.0: + updated["eta"] = _format_float(p_capture / float(p_opt[idx])) + else: + updated["eta"] = "" + else: + updated["eta"] = "" + updated_rows.append(updated) + return updated_rows + + +def _write_csv_exact(csv_path: Path, fieldnames: list[str], rows: list[dict[str, str]]) -> None: + with csv_path.open("w", newline="") as fh: + writer = csv.DictWriter(fh, fieldnames=fieldnames) + writer.writeheader() + writer.writerows(rows) + + +def _masked_summary(before_rows: list[dict[str, str]], after_rows: list[dict[str, str]]) -> tuple[int, int, int]: + before = [str(row.get("masked", "false")).strip().lower() == "true" for row in before_rows] + after = [str(row.get("masked", "false")).strip().lower() == "true" for row in after_rows] + changed = sum(1 for b, a in zip(before, after) if b != a) + return sum(before), sum(after), changed + + +def retabulate_csv(csv_path: Path, h5_path: Path, *, write: bool) -> tuple[int, int, int]: + fieldnames, rows = _load_csv_exact(csv_path) + missing = [column for column in REQUIRED_COLUMNS if column not in fieldnames] + if missing: + raise RuntimeError(f"CSV missing required columns {missing}: {csv_path}") + + updated_rows = _retabulate_rows(rows, h5_path) + summary = _masked_summary(rows, updated_rows) + if write: + _write_csv_exact(csv_path, fieldnames, updated_rows) + return summary + + +def iter_targets(repo: Path) -> Iterator[tuple[str, int, Path, Path]]: + """Yield (label, flap angle, CSV path, hinge-H5 path) retabulation targets.""" + for angle, meta in sorted(FLAPS.items()): + h5_path = repo / meta["h5"] + for label, rel_dir in CSV_GROUPS: + csv_path = repo / rel_dir / f"capture_efficiency_VGM{angle}.csv" + yield label, angle, csv_path, h5_path + + +def parse_args() -> argparse.Namespace: + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument( + "--repo", + default=str(Path(__file__).resolve().parents[1]), + help="Repository root (default: parent of scripts/)", + ) + parser.add_argument( + "--dry-run", + action="store_true", + help="Print masked-change summary without writing CSVs", + ) + return parser.parse_args() + + +def main() -> int: + args = parse_args() + repo = Path(args.repo).resolve() + saw_target = False + + for label, angle, csv_path, h5_path in iter_targets(repo): + if not csv_path.exists(): + print(f"[skip] {label} VGM-{angle}: missing CSV {csv_path}") + continue + if not h5_path.exists(): + print(f"[skip] {label} VGM-{angle}: missing H5 {h5_path}") + continue + saw_target = True + masked_before, masked_after, changed = retabulate_csv(csv_path, h5_path, write=not args.dry_run) + verb = "would update" if args.dry_run else "updated" + print( + f"[{verb}] {label} VGM-{angle}: masked {masked_before} -> {masked_after} " + f"({changed} rows changed)" + ) + + if not saw_target: + print("ERROR: No target CSV/H5 pairs found") + return 2 + return 0 + + +if __name__ == "__main__": + raise SystemExit(main()) diff --git a/tests/test_hydro_retabulation.py b/tests/test_hydro_retabulation.py new file mode 100644 index 0000000..0321db8 --- /dev/null +++ b/tests/test_hydro_retabulation.py @@ -0,0 +1,74 @@ +import csv +import sys +import tempfile +import unittest +from pathlib import Path + +import numpy as np + +REPO_ROOT = Path(__file__).resolve().parents[1] +SCRIPTS_DIR = REPO_ROOT / "scripts" +if str(SCRIPTS_DIR) not in sys.path: + sys.path.insert(0, str(SCRIPTS_DIR)) + +import passive_vs_optpassive_sweep # noqa: E402 +import retabulate_hydro_columns # noqa: E402 + + +class HydroRetabulationTests(unittest.TestCase): + def test_popt_curve_from_h5_clamps_b55_non_negative_for_all_hinged_flaps(self) -> None: + for angle, meta in passive_vs_optpassive_sweep.FLAPS.items(): + with self.subTest(angle=angle): + _, b55, _, _, _ = passive_vs_optpassive_sweep.popt_curve_from_h5( + REPO_ROOT / meta["h5"], + passive_vs_optpassive_sweep.PERIOD_GRID, + ) + self.assertGreaterEqual(float(np.min(b55)), 0.0) + + def test_hinge_basis_mask_differs_from_cg_basis_for_vgm45(self) -> None: + periods_s = np.array([1.75, 3.5], dtype=float) + _, b55_hinge, _, _, masked_hinge = passive_vs_optpassive_sweep.popt_curve_from_h5( + REPO_ROOT / passive_vs_optpassive_sweep.FLAPS[45]["h5"], + periods_s, + ) + _, b55_cg, _, _, masked_cg = passive_vs_optpassive_sweep.popt_curve_from_h5( + REPO_ROOT / "hydroData" / "vgoswec_45.h5", + periods_s, + ) + + self.assertFalse(bool(masked_hinge[0])) + self.assertTrue(bool(masked_cg[0])) + self.assertNotEqual(bool(masked_hinge[0]), bool(masked_cg[0])) + self.assertGreater(float(b55_hinge[0]), passive_vs_optpassive_sweep.MASK_B55_THRESHOLD) + self.assertLessEqual(float(b55_cg[0]), passive_vs_optpassive_sweep.MASK_B55_THRESHOLD) + self.assertGreater(float(b55_hinge[1]), float(b55_cg[1])) + + def test_retabulation_preserves_p_capture_column_text(self) -> None: + source_csv = REPO_ROOT / "analysis" / "passive" / "capture_efficiency_VGM45.csv" + h5_path = REPO_ROOT / passive_vs_optpassive_sweep.FLAPS[45]["h5"] + + with tempfile.TemporaryDirectory() as tmpdir: + csv_copy = Path(tmpdir) / source_csv.name + csv_copy.write_text(source_csv.read_text()) + + with csv_copy.open(newline="") as fh: + before_reader = csv.DictReader(fh) + before_fieldnames = list(before_reader.fieldnames or []) + before_rows = list(before_reader) + + retabulate_hydro_columns.retabulate_csv(csv_copy, h5_path, write=True) + + with csv_copy.open(newline="") as fh: + after_reader = csv.DictReader(fh) + after_fieldnames = list(after_reader.fieldnames or []) + after_rows = list(after_reader) + + self.assertEqual(before_fieldnames, after_fieldnames) + self.assertEqual( + [row["P_capture_W"] for row in before_rows], + [row["P_capture_W"] for row in after_rows], + ) + + +if __name__ == "__main__": + unittest.main()