Mathematical analysis of the pressure-solver discontinuity (issue #77)

discontinuity_analysis/.gitignore

@@ -0,0 +1,8 @@
+# latexmk / TeX intermediates
+*.aux
+*.log
+*.out
+*.toc
+*.fls
+*.fdb_latexmk
+*.synctex.gz

discontinuity_analysis/README.md

@@ -0,0 +1,49 @@
+# Discontinuity analysis of the tcpyPI pressure solver (issue #77 / PR #77)
+
+**The definitive writeup.** `discontinuity_analysis.pdf` (build:
+`latexmk -pdf`) is a self-contained mathematical account of the period-2
+convergence failure in the `pcmin` fixed-point iteration, written for a
+reader with no atmospheric-physics background. Every headline number was
+additionally reproduced by a second, independently written implementation of
+the map and CAPE routines.
+
+Summary of findings:
+
+- The failing sounding is the environment of **Hurricane Francine** (18 UTC
+  2024-09-11, Gulf of Mexico, hours before landfall; NHC best track confirms
+  the position) under a −1.6 K counterfactual SST.
+- The fixed-point map `g` is discontinuous; its only diagonal crossing falls
+  inside a 0.837 hPa jump at x_d ≈ 950.8289947515 hPa, so **no fixed point
+  exists** and the iteration locks into a superattracting, bit-exact 2-cycle.
+- The jump comes from one **discontinuous selection rule** in `cape()`
+  (integrate signed buoyancy to the topmost positive grid level, clamp at
+  zero): a one-bit sign decision at 350 hPa — on a feature worth
+  ~1e-2 J/kg, and ballistically unreachable (W = −61 J/kg there) — toggles
+  287 J/kg of negative area and flips eyewall CAPE between 0 and 224 J/kg.
+  Not quadrature, not grid resolution, not interpolation smoothness (each
+  tested and refuted), not floating point.
+- Smooth toy models separate the genuine discontinuities of the underlying
+  variational problem (argmax jumps with continuous value — Berge; reachable
+  energy at barrier tangencies, which is free-convection physics, not PI's
+  forced eyewall ascent) from the code's artificial one.
+- Jumps of `g` are generic (20% of 1,200 sampled real columns); the fatal
+  diagonal-in-gap placement is rare (~0.04 K of SST for this column) but is
+  found systematically by exactly the attribution-style parameter sweeps that
+  produced the case.
+- Repair: CAPE = max_q W(q) (value function; same quadrature and candidates,
+  one changed selection) + residual-checked bracketed root-finding with the
+  explicit non-hypercane endpoint hypothesis. Exact corrected roots:
+  x* = 950.3136889226 hPa (VMAX 83.9065 m/s, PMIN 920.7517 hPa) for #77 and
+  x* = 960.7687 hPa (VMAX 74.4149 m/s) for the 1995 case. On the
+  7,093-column benchmark: zero flag changes, 98.1% of VMAX within 0.1 m/s,
+  1.9% fragile-tail changes up to 7.9 m/s. The period-2 "rescue" becomes
+  dead code (its midpoint was never a root).
+
+Layout:
+
+- `discontinuity_analysis.tex` / `.pdf` — the document (build: `latexmk -pdf`)
+- `figures/` — all figures (PDF used by the doc, PNG previews)
+- `scripts/` — the analysis scripts behind every number and figure; data
+  files (`*.npz`) are regenerated by running `datagen.py`, `humphunt.py`,
+  `refine.py`, `fixdata.py`, `validate.py` against branch
+  `upstream/modernize_eg`

discontinuity_analysis/scripts/caseC2.py

@@ -0,0 +1,55 @@
+"""Regenerate the 'healthy column' case with a representative interior-LNB
+column (m=0, j=20, i=1). Saves healthy2.npz."""
+
+import os
+import sys
+
+SCRATCH = os.path.dirname(os.path.abspath(__file__))
+sys.path.insert(0, os.path.join(SCRATCH, "branch", "src"))
+sys.path.insert(0, SCRATCH)
+
+import numpy as np
+import xarray as xr
+from tcpyPI import utilities
+
+from datagen import g_scan, iterate_raw, W_curve
+
+M, J, I = 0, 20, 1
+
+ds = xr.open_dataset(os.path.join(SCRATCH, "sample_data.nc"))
+p_s = ds["p"].values.astype(float)
+sst = float(ds["sst"].values[M, J, I])
+msl = float(ds["msl"].values[M, J, I])
+T = utilities.T_Ctok(ds["t"].values[M, :, J, I].astype(float))
+R = ds["r"].values[M, :, J, I].astype(float) * 0.001
+R[np.isnan(R)] = 0.0
+lat = float(ds["lat"].values[J])
+lon = float(ds["lon"].values[I])
+month = int(ds["month"].values[M]) if "month" in ds else M + 1
+print(f"column: month={month} lat={lat} lon={lon} SST={sst:.2f}")
+
+SSTK = utilities.T_Ctok(sst)
+ES0 = utilities.es_cc(sst)
+
+save = {}
+pms = np.arange(850.0, 1005.0, 0.02)
+sc = g_scan(pms, SSTK, msl, p_s, T, R, ES0)
+save["c_pm"] = pms
+save["c_g"] = sc[:, 0]
+st, xs = iterate_raw(SSTK, msl, p_s, T, R, ES0, nmax=60)
+save["c_iter"] = xs
+print(f"iteration: status={st}, fixed point ~ {xs[-1]:.4f}, n={len(xs)}")
+
+pmc = xs[-1]
+PP = min(pmc, 1000.0)
+pts, Ws, b, Pt = W_curve(SSTK, utilities.rv(ES0, PP), PP, T, R, p_s)
+save["c_W_p"] = pts
+save["c_W_W"] = Ws
+save["c_b"] = b
+save["c_blevels"] = p_s[: int(np.count_nonzero(p_s > 50))]
+kmax = int(np.argmax(Ws))
+print(f"argmax LNB: {pts[kmax]:.2f} hPa, CAPE={Ws[kmax]:.1f}; "
+      f"b sign change: b[{save['c_blevels'][-2]:.0f}]={b[-2]:+.2f}, "
+      f"b[{save['c_blevels'][-1]:.0f}]={b[-1]:+.2f}")
+np.savez(os.path.join(SCRATCH, "healthy2.npz"), **save)
+print("saved healthy2.npz")

discontinuity_analysis/scripts/compare.py

@@ -0,0 +1,52 @@
+"""Compare branch (pcmin CAPE + rescue) vs branch_fix (max-W CAPE, no rescue)."""
+
+import numpy as np
+import xarray as xr
+
+a = np.load("branch_pi.npz")      # baseline: branch as merged
+b = np.load("branch_fix_pi.npz")  # max-W CAPE, rescue deleted
+
+both1 = (a["IFL"] == 1) & (b["IFL"] == 1)
+n_both = int(both1.sum())
+
+# flag agreement
+flag_diff = np.argwhere(a["IFL"] != b["IFL"])
+print(f"columns with IFL=1 in both: {n_both}")
+print(f"columns where IFL differs: {len(flag_diff)}")
+for idx in flag_diff[:20]:
+    m, j, i = idx
+    print(f"  (m={m},j={j},i={i}): IFL {a['IFL'][m,j,i]} -> {b['IFL'][m,j,i]}, "
+          f"VMAX {a['VMAX'][m,j,i]:.4f} -> {b['VMAX'][m,j,i]:.4f}")
+
+for var in ["VMAX", "PMIN", "TO", "OTL"]:
+    x, y = a[var], b[var]
+    d = np.abs(x - y)[both1]
+    ident = int((d == 0).sum())
+    print(f"{var}: bit-identical {ident}/{n_both} "
+          f"({100*ident/n_both:.2f}%), max|d|={np.nanmax(d):.6g}, "
+          f"n(|d|>1e-8)={int((d>1e-8).sum())}, n(|d|>0.1)={int((d>0.1).sum())}")
+
+# which columns changed materially, and were they rescue cases?
+d_vmax = np.abs(a["VMAX"] - b["VMAX"])
+chg = np.argwhere(both1 & (d_vmax > 1e-8))
+print(f"\ncolumns with VMAX change > 1e-8: {len(chg)}")
+for idx in chg[:30]:
+    m, j, i = idx
+    print(f"  (m={m},j={j},i={i}): VMAX {a['VMAX'][m,j,i]:.6f}->{b['VMAX'][m,j,i]:.6f} "
+          f"TO {a['TO'][m,j,i]:.3f}->{b['TO'][m,j,i]:.3f} "
+          f"OTL {a['OTL'][m,j,i]:.2f}->{b['OTL'][m,j,i]:.2f}")
+
+# MATLAB reference comparison for both
+ml = xr.open_dataset("matlab_ref.nc")
+print("\nMATLAB ref vars:", list(ml.data_vars))
+mv = ml["Vmax"].values
+mp = ml["Pmin"].values
+mfl = ml["PI_flag"].values if "PI_flag" in ml else None
+print("MATLAB shapes:", mv.shape)
+
+for name, r in [("branch", a), ("fix", b)]:
+    ok = (r["IFL"] == 1) & np.isfinite(mv) & (mfl == 1 if mfl is not None else True)
+    dv = np.abs(r["VMAX"] - mv)[ok]
+    dp = np.abs(r["PMIN"] - mp)[ok]
+    print(f"{name}: vs MATLAB over {int(ok.sum())} cols: "
+          f"VMAX max|d|={np.nanmax(dv):.6f} PMIN max|d|={np.nanmax(dp):.6f}")