make_figures.py

"""
make_figures.py — generate the validation figures for the README.

Produces three PNGs in ./figures/ :
  1. operating_point.png  — pump curve vs system curve, operating point marked
  2. transient.png        — BDF transient approaching the analytical steady state
  3. validation_rth.png   — model thermal resistance vs the closed-form reference

Run after loop1d.py is importable from the same directory:
    python make_figures.py
"""

import os
import numpy as np
import matplotlib.pyplot as plt

from loop1d import (Loop, solve_hydraulic, solve_transient,
                    solve_steady, steady_check, hydraulic_check)

os.makedirs("figures", exist_ok=True)
plt.rcParams.update({
    "figure.dpi": 130, "savefig.dpi": 130, "font.size": 11,
    "axes.grid": True, "grid.alpha": 0.3, "axes.axisbelow": True,
})

# --------------------------------------------------------------------------
# Figure 1 — hydraulic operating point (pump curve meets system curve)
# --------------------------------------------------------------------------
loop = Loop()
mdot_op = solve_hydraulic(loop)
m = np.linspace(0.0, 0.30, 300)

Ksys = loop.cold_plate.K + loop.hx.K + loop.pipe.K
dp_sys = Ksys * m**2

fig, ax = plt.subplots(figsize=(6.4, 4.2))
for spd, c in [(1.0, "#185FA5"), (0.8, "#1D9E75"), (0.6, "#D85A30")]:
    dp_pump = loop.pump.dp0 * spd**2 - loop.pump.k * m**2
    ax.plot(m, dp_pump / 1e3, color=c, lw=2, label=f"pump curve, speed = {spd:.1f}")
    # operating point at this speed
    mo = np.sqrt(loop.pump.dp0 * spd**2 / (Ksys + loop.pump.k))
    ax.plot(mo, (Ksys * mo**2) / 1e3, "o", color=c, ms=7,
            markeredgecolor="white", zorder=5)

ax.plot(m, dp_sys / 1e3, "k--", lw=1.8, label="system curve (losses)")
ax.set_xlabel("mass flow rate  $\\dot{m}$  [kg/s]")
ax.set_ylabel("pressure rise / drop  $\\Delta p$  [kPa]")
ax.set_title("Hydraulic operating point: pump curve $\\cap$ system curve")
ax.set_ylim(0, 240)
ax.set_xlim(0, 0.30)
ax.legend(frameon=False, fontsize=9.5)
fig.tight_layout()
fig.savefig("figures/operating_point.png")
plt.close(fig)

# --------------------------------------------------------------------------
# Figure 2 — transient approach to the analytical steady state
# --------------------------------------------------------------------------
loop = Loop()
mdot = solve_hydraulic(loop)
sol = solve_transient(loop, mdot, t_end=60.0)
Tw_ss, Tcp_ss, Thx_ss, _ = steady_check(loop, mdot)

t = np.linspace(0, 60, 400)
Y = sol.sol(t)

fig, ax = plt.subplots(figsize=(6.4, 4.2))
labels = ["wall (junction proxy)", "cold-plate fluid", "HX-outlet fluid"]
colors = ["#993C1D", "#185FA5", "#1D9E75"]
refs = [Tw_ss, Tcp_ss, Thx_ss]
for i in range(3):
    ax.plot(t, Y[i], color=colors[i], lw=2, label=labels[i])
    ax.axhline(refs[i], color=colors[i], ls=":", lw=1.2, alpha=0.8)
ax.set_xlabel("time  [s]")
ax.set_ylabel("temperature  [°C]")
ax.set_title("Transient (BDF) settling onto the analytical steady state")
ax.text(0.97, 0.5, "dotted lines = closed-form steady values",
        transform=ax.transAxes, ha="right", fontsize=9, color="#444441")
ax.set_xlim(0, 60)
ax.legend(frameon=False, fontsize=9.5, loc="center right",
          bbox_to_anchor=(1.0, 0.32))
fig.tight_layout()
fig.savefig("figures/transient.png")
plt.close(fig)

# --------------------------------------------------------------------------
# Figure 3 — thermal-resistance validation vs closed form, across flow
# --------------------------------------------------------------------------
speeds = np.linspace(0.45, 1.0, 12)
mdots, R_model, R_ref = [], [], []
for spd in speeds:
    lp = Loop()
    lp.pump.speed = spd
    md = solve_hydraulic(lp)
    Tw, Tcp, Thx = solve_steady(lp, md)
    mdots.append(md)
    R_model.append((Tw - Thx) / lp.Q_chip)                 # from the solver
    R_ref.append(1.0 / lp.cold_plate.hA + 1.0 / (md * lp.fluid.cp))  # closed form

mdots = np.array(mdots)
fig, ax = plt.subplots(figsize=(6.4, 4.2))
ax.plot(mdots, np.array(R_ref) * 1e3, "-", color="#185FA5", lw=2,
        label=r"closed form  $1/hA + 1/(\dot{m}c_p)$")
ax.plot(mdots, np.array(R_model) * 1e3, "o", color="#993C1D", ms=6,
        markeredgecolor="white", label="model (steady solve)")
ax.set_xlabel("mass flow rate  $\\dot{m}$  [kg/s]")
ax.set_ylabel("cold-plate thermal resistance  $R_{th}$  [mK/W]")
ax.set_title("Validation: solver output vs analytical reference")
ax.legend(frameon=False, fontsize=10)
fig.tight_layout()
fig.savefig("figures/validation_rth.png")
plt.close(fig)

# residual report for the README badge / text
err = np.max(np.abs(np.array(R_model) - np.array(R_ref)) / np.array(R_ref))
print(f"Figures written to ./figures/")
print(f"Max R_th relative error (model vs closed form): {err:.2e}")