Add config-driven reproducibility layer

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

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

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()

requirements.txt

@@ -1,4 +1,5 @@
 numpy>=1.24
 scipy>=1.10
 matplotlib>=3.7
+pyyaml>=6.0
 pytest>=7.0

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
+  }
+}