plotConvergence.py

#!/usr/bin/env python3
"""
Grid-convergence post-processing for the K&O dam-break-with-obstacle study.

Reads study_coarse / medium / fine / veryfine and reports, for each mesh:
  - surge-front curve (Z/L vs T*)         [kinematic]
  - peak obstacle force F_x               [converges]
  - peak impact pressure                  [diverges - local singularity]
  - impact IMPULSE  J = integral p dt      [bounded -> converges]
    over a FIXED window (same bounds for all meshes, so it is a fair test).

Figures:  front_overlay.png ,  convergence.png (force | peak P | impulse)
Run from the study root:  python3 plotConvergence.py     (no re-running needed)
"""
import os
import glob
import math

import numpy as np
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt

L, H, G, RHO = 0.146, 0.292, 9.81, 1000.0
AREA = 4 * L * 4 * L
IMPACT_WINDOW = (0.10, 0.30)     # [s] brackets the first impact on every mesh
HERE = os.path.dirname(os.path.abspath(__file__))
LEVELS = ["coarse", "medium", "fine", "veryfine"]


def n_cells(case):
    log = os.path.join(case, "log.checkMesh")
    if not os.path.isfile(log):
        return None
    for line in open(log):
        s = line.strip()
        if s.startswith("cells:"):
            return int(s.split()[1])
    return None


def peak_force(case):
    f = glob.glob(os.path.join(case, "postProcessing", "obstacleForces", "*", "force*.dat"))
    if not f:
        return None
    fx = []
    for line in open(sorted(f)[0]):
        if line.startswith("#") or not line.strip():
            continue
        toks = line.replace("(", " ").replace(")", " ").split()
        try:
            fx.append(float(toks[1]))
        except (ValueError, IndexError):
            pass
    return max(fx) if fx else None


def _probe_series(case):
    """Return (t, p_lowest_probe) for the y=4 mm probe, or None."""
    f = glob.glob(os.path.join(case, "postProcessing", "pressureProbes", "*", "p"))
    if not f:
        return None
    d = np.loadtxt(sorted(f)[0], comments="#")
    return d[:, 0], d[:, 1]


def peak_pressure(case):
    s = _probe_series(case)
    return s[1].max() if s else None


def impulse(case, window=IMPACT_WINDOW):
    """Time-integral of pressure over the fixed impact window [Pa.s]."""
    s = _probe_series(case)
    if not s:
        return None
    t, p = s
    m = (t >= window[0]) & (t <= window[1])
    if m.sum() < 2:
        return None
    _trap = getattr(np, "trapezoid", getattr(np, "trapz", None))
    return float(_trap(p[m], t[m]))


def front_curve(case):
    files = glob.glob(os.path.join(case, "postProcessing", "frontTracking",
                                   "*", "floorLine_*.xy"))
    pts = []
    for fp in files:
        t = float(os.path.basename(os.path.dirname(fp)))
        a = np.loadtxt(fp)
        if a.ndim == 1:
            a = a[None, :]
        wet = np.where(a[:, 1] >= 0.5)[0]
        xf = a[wet, 0].max() if wet.size else 0.0
        pts.append((t, xf))
    if not pts:
        return None
    pts = np.array(sorted(pts))
    return pts[:, 0] * math.sqrt(2 * G / L), pts[:, 1] / L


# ---- gather ----
rows = []
for lvl in LEVELS:
    case = os.path.join(HERE, f"study_{lvl}")
    rows.append(dict(level=lvl, case=case, N=n_cells(case),
                     Fx=peak_force(case), P=peak_pressure(case), J=impulse(case)))
rows = [r for r in rows if r["N"]]
if not rows:
    raise SystemExit("No study_* results found - run the cases first.")
for r in rows:
    r["h"] = math.sqrt(AREA / r["N"])
rows.sort(key=lambda x: x["N"])

# ---- figure 1: front overlay ----
plt.figure(figsize=(7, 5))
for r in rows:
    fc = front_curve(r["case"])
    if fc:
        plt.plot(fc[0], fc[1], lw=2, label=f"{r['level']} (N={r['N']})")
plt.axhline(4.0, color="0.7", ls="--", lw=1)
plt.xlabel(r"$T^* = t\sqrt{2g/L}$"); plt.ylabel(r"$Z/L$")
plt.title("Surge front - mesh sensitivity")
plt.grid(alpha=0.3); plt.legend()
plt.tight_layout(); plt.savefig(os.path.join(HERE, "front_overlay.png"), dpi=150)

# ---- figure 2: force | peak P | impulse vs h ----
h  = np.array([r["h"] for r in rows])
Fx = np.array([r["Fx"] for r in rows])
P  = np.array([r["P"] for r in rows]) / 1000.0
J  = np.array([r["J"] for r in rows])

fig, ax = plt.subplots(1, 3, figsize=(15, 4.5))
ax[0].plot(h, Fx, "o-", lw=2);            ax[0].set_title("Peak force (converges)")
ax[0].set_ylabel(r"peak $F_x$ [N]")
ax[1].plot(h, P, "s-", lw=2, color="C3"); ax[1].set_title("Peak pressure (diverges)")
ax[1].set_ylabel("peak P [kPa]")
ax[2].plot(h, J, "^-", lw=2, color="C2"); ax[2].set_title("Impact impulse (converges)")
ax[2].set_ylabel(r"$\int p\,dt$ over impact [Pa$\cdot$s]")
for a in ax:
    a.set_xlabel("cell size h [m]"); a.grid(alpha=0.3); a.invert_xaxis()
plt.tight_layout(); plt.savefig(os.path.join(HERE, "convergence.png"), dpi=150)

# ---- table ----
print(f"impact window for impulse: {IMPACT_WINDOW[0]}-{IMPACT_WINDOW[1]} s\n")
print(f"{'level':9s}{'N cells':>9s}{'h[m]':>9s}"
      f"{'peakFx[N]':>11s}{'peakP[kPa]':>12s}{'impulse[Pa.s]':>15s}")
prev = None
for r in rows:
    line = (f"{r['level']:9s}{r['N']:>9d}{r['h']:>9.4f}"
            f"{r['Fx']:>11.3f}{r['P']/1000:>12.3f}{r['J']:>15.3f}")
    if prev:
        dF = 100*abs(r["Fx"]-prev["Fx"])/abs(r["Fx"])
        dP = 100*abs(r["P"] -prev["P"]) /abs(r["P"])
        dJ = 100*abs(r["J"] -prev["J"]) /abs(r["J"])
        line += f"  (dFx {dF:.1f}% dP {dP:.1f}% dJ {dJ:.1f}%)"
    print(line)
    prev = r
print("\nPeak pressure %-change should keep growing; impulse %-change should "
      "shrink. That contrast is the headline result.")
print("Wrote front_overlay.png and convergence.png")