"""Threshold calibration for the mixed-derivative ADI benchmark. Usage: python scheme_survey.py This is the evidence behind the acceptance thresholds in grade.py. The task permits ANY ADI scheme that is stable and second order in the presence of the mixed derivative, so the rubric must not be tuned to one particular method. We run the whole in 't Hout-Foulon ADI family on the benchmark problem -- Douglas (Do), Craig-Sneyd (CS), Modified Craig-Sneyd (MCS) and Hundsdorfer-Verwer (HV), at several values of the parameter theta -- together with the classic failure modes, and check each one against the thresholds grade.py actually enforces. Expected outcome: every row expected to be second order on this mesh sequence is accepted, every other row rejected. The script exits non-zero on any mismatch, so the coverage claim below cannot drift away from what grade.py actually enforces. Background (in 't Hout & Welfert, 2007/2009), with a mixed derivative in 2D: * Do is only FIRST order, for any theta; * CS is second order ONLY at theta = 1/2 -- away from it the measured order falls to ~1.2 (theta = 0.8 -> 1.205, theta = 2 -> 1.410), and at theta = 1/3 it is unstable and overflows; * MCS is second order and unconditionally stable for theta >= 1/3; * HV is second order and unconditionally stable for theta >= 1/2 + sqrt(3)/6. Note the distinction this script exists to enforce: MCS and HV are second order *asymptotically* for any admissible theta, but at large theta the asymptotic regime lies beyond h = 0.00625, so the order OBSERVED on the prescribed meshes stops being ~2. The band ends at MCS theta = 10 (79.9x total reduction, fitted order 1.652 -- short of the published 100x and [1.70, 2.30]). MCS at theta = 8 IS accepted: its pointwise orders are an untidy transient, [0.62, 1.82, 2.45, 1.92], but its reduction (112x) and fitted order (1.788) clear the published bar. The task requires the observed order to approach 2 on exactly h = 0.1*2^-p, p = 0..4, so the accepted band is MEASURED here, not inherited from the stability theory -- which is why this script fails loudly rather than documenting a claim that could drift. """ import os import sys import numpy as np from scipy.linalg import solve_banded sys.path.insert(0, os.path.dirname(os.path.abspath(__file__))) from grade import ( # noqa: E402 (thresholds are the point of this script) FINEST_MAX_ERR, LS_ORDER_HI, LS_ORDER_LO, MIN_TOTAL_REDUCTION, PT_ORDER_HI, PT_ORDER_LO, TIER1_MAX_ERR, ) A_COEF, B_COEF, C_COEF = 1.0, 0.8, 0.5 T_FINAL = 1.0 MODES = [(1.0, 0.7, 1.3, 0.0), (0.7, 1.4, -0.6, 0.3), (0.5, 1.0, 0.8, 0.7)] H_VALUES = [0.1 * 2.0**(-p) for p in range(5)] def u_exact(x, y, t): s = np.zeros_like(x * y * 1.0) for (amp, al, be, ph) in MODES: lam = A_COEF * al * al + 2 * C_COEF * al * be + B_COEF * be * be s = s + amp * np.exp(-lam * t) * np.sin(al * x + be * y + ph) return s def run(h, scheme, theta, T=T_FINAL): """One ADI step family, parameterised by scheme name and theta.""" m = int(round(1.0 / h)) k = h N = int(round(T / k)) h2 = h * h ax = np.linspace(0.0, 1.0, m + 1) X, Y = np.meshgrid(ax, ax, indexing="ij") def fill(V, t): G = u_exact(X, Y, t).copy() G[1:m, 1:m] = V return G def Dxx(G): return (G[2:m + 1, 1:m] - 2 * G[1:m, 1:m] + G[0:m - 1, 1:m]) / h2 def Dyy(G): return (G[1:m, 2:m + 1] - 2 * G[1:m, 1:m] + G[1:m, 0:m - 1]) / h2 def Dxy(G): return (G[2:m + 1, 2:m + 1] - G[2:m + 1, 0:m - 1] - G[0:m - 1, 2:m + 1] + G[0:m - 1, 0:m - 1]) / (4 * h2) def F0(V, t): return 2 * C_COEF * Dxy(fill(V, t)) def F1(V, t): return A_COEF * Dxx(fill(V, t)) def F2(V, t): return B_COEF * Dyy(fill(V, t)) def F(V, t): G = fill(V, t) return 2 * C_COEF * Dxy(G) + A_COEF * Dxx(G) + B_COEF * Dyy(G) rx = theta * k * A_COEF / h2 ry = theta * k * B_COEF / h2 abx = np.zeros((3, m - 1)); abx[0, 1:] = -rx; abx[1, :] = 1 + 2 * rx; abx[2, :-1] = -rx aby = np.zeros((3, m - 1)); aby[0, 1:] = -ry; aby[1, :] = 1 + 2 * ry; aby[2, :-1] = -ry def solve_x(R, t_imp): gb = u_exact(X, Y, t_imp) Rb = R.copy() Rb[0, :] += rx * gb[0, 1:m] Rb[-1, :] += rx * gb[m, 1:m] return solve_banded((1, 1), abx, Rb) def solve_y(R, t_imp): gb = u_exact(X, Y, t_imp) Rb = R.copy() Rb[:, 0] += ry * gb[1:m, 0] Rb[:, -1] += ry * gb[1:m, m] return solve_banded((1, 1), aby, Rb.T).T U = u_exact(X, Y, 0.0)[1:m, 1:m].copy() for n in range(N): tn, tn1 = n * k, (n + 1) * k F0n, F1n, F2n = F0(U, tn), F1(U, tn), F2(U, tn) Fn = (F1n + F2n) if scheme == "drop" else (F0n + F1n + F2n) Y0 = U + k * Fn # explicit predictor Y1 = solve_x(Y0 - theta * k * F1n, tn1) # x sweep Y2 = solve_y(Y1 - theta * k * F2n, tn1) # y sweep if scheme in ("Do", "drop"): U = Y2 continue if scheme == "CS": Yh0 = Y0 + 0.5 * k * (F0(Y2, tn1) - F0n) elif scheme == "MCS": Yh0 = (Y0 + theta * k * (F0(Y2, tn1) - F0n) + (0.5 - theta) * k * (F(Y2, tn1) - Fn)) elif scheme == "HV": Yh0 = Y0 + 0.5 * k * (F(Y2, tn1) - Fn) else: raise ValueError(scheme) if scheme == "HV": # HV subtracts F_j evaluated at (Y2, t_{n+1}) Yh1 = solve_x(Yh0 - theta * k * F1(Y2, tn1), tn1) U = solve_y(Yh1 - theta * k * F2(Y2, tn1), tn1) else: # CS / MCS subtract F_j at (U, t_n) Yh1 = solve_x(Yh0 - theta * k * F1n, tn1) U = solve_y(Yh1 - theta * k * F2n, tn1) Ufull = u_exact(X, Y, T).copy() Ufull[1:m, 1:m] = U return float(np.max(np.abs(Ufull - u_exact(X, Y, T)))) CASES = [ # (label, scheme, theta, expected to be accepted?) # --- second order on this mesh sequence: MUST be accepted --------------- ("MCS theta=1/3", "MCS", 1.0 / 3.0, True), ("MCS theta=1/2", "MCS", 0.5, True), ("MCS theta=0.8", "MCS", 0.8, True), ("MCS theta=1", "MCS", 1.0, True), ("MCS theta=1.5", "MCS", 1.5, True), ("MCS theta=2", "MCS", 2.0, True), ("MCS theta=3", "MCS", 3.0, True), ("MCS theta=5", "MCS", 5.0, True), ("MCS theta=8", "MCS", 8.0, True), ("CS theta=1/2", "CS", 0.5, True), ("HV theta=1/2", "HV", 0.5, True), ("HV theta=1/2+sqrt(3)/6", "HV", 0.5 + np.sqrt(3) / 6, True), ("HV theta=1", "HV", 1.0, True), ("HV theta=2", "HV", 2.0, True), ("HV theta=3", "HV", 3.0, True), ("HV theta=5", "HV", 5.0, True), ("HV theta=8", "HV", 8.0, True), # --- NOT second order on this mesh sequence: must be rejected ----------- # MCS/HV are second order asymptotically for any theta, but at large theta # the asymptotic regime lies beyond h = 0.00625, so the observed order on the # PRESCRIBED meshes is not ~2: MCS at theta = 10 achieves only a 79.9x total # reduction and a fitted order of 1.652 over the five meshes. The instruction # requires the order to approach 2 on exactly this sequence, so these are # correctly out. This is where the accepted band genuinely ends. ("MCS theta=10 (pre-asymptotic)", "MCS", 10.0, False), ("MCS theta=20 (pre-asymptotic)", "MCS", 20.0, False), # --- genuine failure modes --------------------------------------------- ("CS theta=1/3 (unstable)", "CS", 1.0 / 3.0, False), ("CS theta=0.8 (1st order)", "CS", 0.8, False), ("CS theta=2 (1st order)", "CS", 2.0, False), ("Do theta=1/2 (1st order)", "Do", 0.5, False), ("Do theta=1 (1st order)", "Do", 1.0, False), ("mixed term dropped", "drop", 0.5, False), ] def verdict(errs): """Apply exactly the criteria grade.py enforces.""" e = np.asarray(errs, float) if not np.all(np.isfinite(e)) or not np.all(e > 0): return False, "non-finite / non-positive" if e[0] >= TIER1_MAX_ERR: return False, "Tier 1 cap" if not np.all(np.diff(e) < 0): return False, "not monotone" if e[0] / e[-1] < MIN_TOTAL_REDUCTION: return False, "total reduction" ls = float(-np.polyfit(np.arange(len(e)), np.log2(e), 1)[0]) if not (LS_ORDER_LO <= ls <= LS_ORDER_HI): return False, "fitted order" pt = np.log2(e[:-1] / e[1:])[-2:] if not all(PT_ORDER_LO <= v <= PT_ORDER_HI for v in pt): return False, "pointwise order" if e[-1] >= FINEST_MAX_ERR: return False, "finest cap" return True, "" def main(): print(f"thresholds: e(0.1) < {TIER1_MAX_ERR:.0e} | e(finest) < {FINEST_MAX_ERR:.0e} | " f"reduction >= {MIN_TOTAL_REDUCTION:.0f}x | fitted order in " f"[{LS_ORDER_LO}, {LS_ORDER_HI}] | last two orders in " f"[{PT_ORDER_LO}, {PT_ORDER_HI}]\n") hdr = (f"{'scheme':<30}{'e(0.1)':>11}{'e(finest)':>12}{'reduction':>10}" f"{'fitted':>8}{'last two orders':>18} {'expected':<13}{'verdict':<8}") print(hdr) print("-" * len(hdr)) ok_all = True for label, sch, th, admissible in CASES: with np.errstate(all="ignore"): errs = [] for h in H_VALUES: try: errs.append(run(h, sch, th)) except Exception: errs.append(float("nan")) e = np.array(errs) accepted, why = verdict(e) red = e[0] / e[-1] if np.all(np.isfinite(e)) and e[-1] > 0 else float("nan") try: ls = float(-np.polyfit(np.arange(5), np.log2(e), 1)[0]) p2 = np.log2(e[:-1] / e[1:])[-2:] except Exception: ls = float("nan"); p2 = np.array([float("nan")] * 2) expected = "2nd-order" if admissible else "not-2nd-ord" got = "accept" if accepted else "reject" agree = (accepted == admissible) ok_all &= agree flag = "" if agree else " <-- MISMATCH" print(f"{label:<30}{e[0]:>11.3e}{e[-1]:>12.3e}{red:>10.1f}{ls:>8.3f}" f"{p2[0]:>9.3f}{p2[1]:>9.3f} {expected:<13}{got:<8}{flag}" + (f" [{why}]" if why else "")) print() if ok_all: print("OK: every scheme that converges at second order on this mesh") print(" sequence is accepted, and everything that does not is rejected.") print() print("COVERAGE CLAIM (this is the precise claim; do not overstate it):") print(" accepted: MCS for theta in [1/3, 8], CS at theta = 1/2,") print(" HV for theta in [1/2, 8].") print(" rejected: schemes that are second order only ASYMPTOTICALLY but") print(" not on h = 0.1*2^-p, p = 0..4 (MCS at theta >= 10),") print(" first-order schemes, and unstable parameter choices.") print(" The instruction requires the observed order to approach 2 on") print(" exactly this mesh sequence, so the second group is correctly out.") else: print("MISMATCH: the thresholds do not separate second-order schemes from") print(" the rest. Recalibrate grade.py.") sys.exit(1) if __name__ == "__main__": main()