diff --git a/README.md b/README.md index 667132e..c2484e3 100644 --- a/README.md +++ b/README.md @@ -66,6 +66,27 @@ print(f"Population served: {med.population_served():,} people") print(f"CO₂ avoided: {med.co2_saved_annual()['total_t_yr']:,} t/yr") ``` +## Reproduce the Published Result (one command) + +All key parameters live in [`config/default.yaml`](config/default.yaml). To +reproduce the figure and metrics in [`results/`](results/) from a fresh clone: + +```bash +git clone https://github.com/defnalk/htgr-desalination.git +cd htgr-desalination +pip install -r requirements.txt +python -m htgr.cli --config config/default.yaml +``` + +Outputs: + +- `results/metrics.json` — reactor transient peaks, thermal-core temperatures, cycle summary, desal production, CO₂ avoided +- `results/htgr_simulation_results.png` — 4-panel summary figure + +To explore variations, copy `config/default.yaml`, edit any parameter +(reactor power, reactivity insertion, MED GOR), and pass the new file via +`--config`. No code changes required. + ## Run the Full Simulation ```bash diff --git a/config/default.yaml b/config/default.yaml new file mode 100644 index 0000000..e12d2ff --- /dev/null +++ b/config/default.yaml @@ -0,0 +1,19 @@ +# Reproducibility config for htgr-desalination. +# Run: python -m htgr.cli --config config/default.yaml + +reactor: + Q_nominal_W: 3.0e8 # 300 MWth nominal thermal power + t_end_s: 400.0 # transient duration + reactivity_pcm: 300 # step reactivity insertion (pcm) + pulse_start_s: 100.0 + pulse_end_s: 150.0 + +desalination: + GOR: 11.1 # gained-output-ratio for the MED train + availability: 0.95 + +output: + results_dir: results + figure_name: htgr_simulation_results.png + metrics_name: metrics.json + random_seed: 0 diff --git a/htgr/cli.py b/htgr/cli.py new file mode 100644 index 0000000..64ebe93 --- /dev/null +++ b/htgr/cli.py @@ -0,0 +1,138 @@ +""" +htgr.cli — config-driven entry point for reproducible HTGR-desal runs. + +Usage: + python -m htgr.cli --config config/default.yaml +""" + +from __future__ import annotations + +import argparse +import json +import random +from pathlib import Path + +import numpy as np +import yaml +from scipy.interpolate import interp1d + +from htgr import PointKinetics, ThermalCore, BraytonRankineCycle, MEDDesalination + + +def load_config(path: Path) -> dict: + with path.open() as f: + return yaml.safe_load(f) + + +def run(config_path: Path) -> dict: + cfg = load_config(config_path) + seed = cfg["output"].get("random_seed", 0) + random.seed(seed) + np.random.seed(seed) + + rcfg = cfg["reactor"] + pk = PointKinetics(Q_nominal=float(rcfg["Q_nominal_W"])) + rho = float(rcfg["reactivity_pcm"]) * 1e-5 + p_start = float(rcfg["pulse_start_s"]) + p_end = float(rcfg["pulse_end_s"]) + t_end = float(rcfg["t_end_s"]) + + def rho_pulse(t): + return rho if p_start <= t <= p_end else 0.0 + + t_kin, P, Q = pk.simulate(t_end=t_end, rho_func=rho_pulse) + + tc = ThermalCore() + Q_interp = interp1d(t_kin, Q, kind="linear", fill_value="extrapolate") + t_th, Tf, Tm, Tc_node = tc.simulate(Q_profile=Q_interp, t_end=t_end) + + cyc = BraytonRankineCycle() + med = MEDDesalination(GOR=float(cfg["desalination"]["GOR"])) + + metrics = { + "config": str(config_path), + "reactor": { + "P_peak_norm": float(P.max()), + "P_peak_t_s": float(t_kin[P.argmax()]), + "Q_peak_MW": float(Q.max() / 1e6), + }, + "thermal": { + "Tf_peak_K": float(Tf.max()), + "Tm_peak_K": float(Tm.max()), + "Tc_peak_K": float(Tc_node.max()), + }, + "cycle_summary": cyc.summary(), + "desal": { + "daily_production_m3": float(med.daily_production()), + "GOR": float(cfg["desalination"]["GOR"]), + }, + "co2_avoided_t_yr": med.co2_saved_annual(), + } + + out_dir = Path(cfg["output"]["results_dir"]) + out_dir.mkdir(parents=True, exist_ok=True) + metrics_path = out_dir / cfg["output"]["metrics_name"] + with metrics_path.open("w") as f: + json.dump(metrics, f, indent=2, default=str) + + _render_figure(cfg, t_kin, P, t_th, Tf, Tm, Tc_node, med, out_dir / cfg["output"]["figure_name"]) + + print(f"✓ metrics → {metrics_path}") + print(f"✓ figure → {out_dir / cfg['output']['figure_name']}") + return metrics + + +def _render_figure(cfg, t_kin, P, t_th, Tf, Tm, Tc_node, med, out_path: Path) -> None: + import matplotlib.pyplot as plt + import matplotlib.gridspec as gridspec + + fig = plt.figure(figsize=(14, 10)) + gs = gridspec.GridSpec(2, 2, hspace=0.40, wspace=0.32) + + ax1 = fig.add_subplot(gs[0, 0]) + ax1.plot(t_kin, P, lw=2.5) + ax1.axvline(cfg["reactor"]["pulse_start_s"], ls=":", color="grey") + ax1.set_xlabel("Time (s)") + ax1.set_ylabel("P / P0") + ax1.set_title("A | Power transient") + ax1.grid(alpha=0.3) + + ax2 = fig.add_subplot(gs[0, 1]) + ax2.plot(t_th, Tf, lw=2, label="Fuel") + ax2.plot(t_th, Tm, lw=2, label="Moderator") + ax2.plot(t_th, Tc_node, lw=2, label="Coolant") + ax2.set_xlabel("Time (s)") + ax2.set_ylabel("T (K)") + ax2.set_title("B | Thermal core") + ax2.legend(fontsize=8) + ax2.grid(alpha=0.3) + + ax3 = fig.add_subplot(gs[1, 0]) + ax3.text(0.05, 0.5, f"Daily production:\n {med.daily_production():.0f} m³/day\nGOR = {cfg['desalination']['GOR']}", + fontsize=14, va="center") + ax3.set_axis_off() + ax3.set_title("C | Desalination") + + ax4 = fig.add_subplot(gs[1, 1]) + co2 = med.co2_saved_annual() + cats = ["Desal", "Elec", "Total"] + vals = [co2["desal_t_yr"] / 1000, co2["elec_t_yr"] / 1000, co2["total_t_yr"] / 1000] + ax4.bar(cats, vals) + ax4.set_ylabel("kt CO₂/yr") + ax4.set_title("D | CO₂ avoided") + ax4.grid(alpha=0.3) + + fig.suptitle("HTGR-desalination — reproducible run", fontweight="bold", y=1.01) + fig.savefig(out_path, dpi=150, bbox_inches="tight") + plt.close(fig) + + +def main() -> None: + parser = argparse.ArgumentParser(prog="python -m htgr.cli") + parser.add_argument("--config", type=Path, default=Path("config/default.yaml")) + args = parser.parse_args() + run(args.config) + + +if __name__ == "__main__": + main() diff --git a/requirements.txt b/requirements.txt index 4944df2..43effb3 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,4 +1,5 @@ numpy>=1.24 scipy>=1.10 matplotlib>=3.7 +pyyaml>=6.0 pytest>=7.0 diff --git a/results/htgr_simulation_results.png b/results/htgr_simulation_results.png new file mode 100644 index 0000000..3592037 Binary files /dev/null and b/results/htgr_simulation_results.png differ diff --git a/results/metrics.json b/results/metrics.json new file mode 100644 index 0000000..73fc78e --- /dev/null +++ b/results/metrics.json @@ -0,0 +1,23 @@ +{ + "config": "config/default.yaml", + "reactor": { + "P_peak_norm": 1.0000000000000078, + "P_peak_t_s": 197.3973973973974, + "Q_peak_MW": 300.00000000000233 + }, + "thermal": { + "Tf_peak_K": 4175.095595719513, + "Tm_peak_K": 2676.0317083271216, + "Tc_peak_K": 1421.9949208643702 + }, + "cycle_summary": "=======================================================\n HTGR Combined Cycle Performance Summary\n=======================================================\n Reactor thermal power : 300.0 MWth\n Helium flow rate : 105.0 kg/s\n Steam flow rate : 12.5 kg/s\n-------------------------------------------------------\n He Brayton net output : -6.0 MWe\n Steam Rankine output : 9.8 MWe\n Combined electrical : 3.8 MWe\n HRSG heat duty : 39.1 MWth\n Cycle thermal efficiency : 1.3%\n-------------------------------------------------------\n Brayton Carnot (950\u2192650\u00b0C): 24.5%\n Rankine Carnot (500\u219270\u00b0C) : 55.6%\n Combined Carnot (950\u219270\u00b0C): 71.9%\n=======================================================", + "desal": { + "daily_production_m3": 11388.6, + "GOR": 11.1 + }, + "co2_avoided_t_yr": { + "desal_t_yr": 20784, + "elec_t_yr": 246331, + "total_t_yr": 267115 + } +} \ No newline at end of file