#!/usr/bin/env python3 """verify_p1p1_tri.py -- L2 verification: triangular P1-P1 Stokes pair. Pure NumPy reference implementation (numpy 2.4.4 at pinning time, 2026-07-10). Pins every number cited in L2 Demo 2, and CORRECTS the prediction recorded in the L1 materials for Exercise 4(c): the pair has EXACT spurious pressure modes on every mesh of the families below; the essential inf-sup constant decays like O(h) after the exact modes are removed. The pair therefore fails through both channels (exact modes, and uniform decay of the essential constant). Setting: Omega = unit square, structured N x N pattern of squares, each split into two triangles; V_h = continuous vector P1 with homogeneous Dirichlet conditions on the whole boundary, Q_h = continuous scalar P1, b(v, q) = -(div v, q). Norms (blueprint v0.6): ||v||_V = ||grad v||_0 (seminorm as norm on (H^1_0)^2), ||q||_0 on Q. Weighted eigenproblem: (B K_V^{-1} B^T) q = lambda M_Q q ; zero eigenvalues <-> ker B^T (exact spurious modes, the constant always among them); beta_tilde := sqrt(first nonzero lambda) is the inf-sup constant of the modified pair (V_h, Q_h ^ complement of the exact modes) -- the 'essential' constant. Pinned facts (asserted below): A. Structured mesh, either diagonal, any N >= 8, any parity: dim ker B^T = 8 (constant + 7 global wave modes). B. Random perturbation of the interior vertices (0.2 h, several seeds), boundary straight: dim ker B^T = 5. The same count holds when the boundary vertices are perturbed as well (generic polygon): 4 of the 7 extra modes are TOPOLOGICAL (survive any geometry with this connectivity), 3 are symmetry artifacts. C. Alternating-diagonal pattern: 8 structured, 7 perturbed -- the topological count depends on the connectivity pattern. D. Essential constant on the structured family decays ~ O(h) (values printed; monotone, ratio N=8 -> N=32 below 0.35). E. Independent matrix-free check: for a computed kernel vector, the patch sums sum_{T ni k} grad(q_h)|_T vanish at interior vertices to ~1e-15 (no assembled matrices involved). F. N = 4 violates the counting inequality (18 velocity dofs < 24 pressure dofs after removing the constant): >= 7 zeros forced; 8 observed. A classification of the topological modes is deliberately left open (L2 exercise sheet, Problem 4(d), open-ended part). """ import numpy as np ME = np.array([[2.0, 1.0, 1.0], [1.0, 2.0, 1.0], [1.0, 1.0, 2.0]]) / 12.0 def build(N, pattern="right", amp_int=0.0, amp_bnd=0.0, seed=0): """Assemble (K_V, B, M_Q) with Dirichlet velocity dofs eliminated. pattern: 'right' (diagonal lower-left -> upper-right), 'left' (Firedrake UnitSquareMesh default), 'alt' (alternating by parity). amp_int / amp_bnd: random vertex perturbation amplitudes, in units of h (boundary flagged by INDEX, so a perturbed boundary yields a generic polygon). """ rng = np.random.default_rng(seed) h = 1.0 / N XY = np.array([[i * h, j * h] for j in range(N + 1) for i in range(N + 1)]) def vid(i, j): return j * (N + 1) + i bnd = np.zeros((N + 1) ** 2, dtype=bool) for j in range(N + 1): for i in range(N + 1): if i in (0, N) or j in (0, N): bnd[vid(i, j)] = True if amp_int > 0: XY[~bnd] += amp_int * h * rng.uniform(-1, 1, size=((~bnd).sum(), 2)) if amp_bnd > 0: XY[bnd] += amp_bnd * h * rng.uniform(-1, 1, size=(bnd.sum(), 2)) cells = [] for j in range(N): for i in range(N): v00, v10 = vid(i, j), vid(i + 1, j) v01, v11 = vid(i, j + 1), vid(i + 1, j + 1) right = (pattern == "right") or (pattern == "alt" and (i + j) % 2 == 0) if right: cells += [(v00, v10, v11), (v00, v11, v01)] else: cells += [(v00, v10, v01), (v10, v11, v01)] nv = (N + 1) ** 2 K = np.zeros((nv, nv)) Mq = np.zeros((nv, nv)) Bx = np.zeros((nv, nv)) By = np.zeros((nv, nv)) for c in cells: (x1, y1), (x2, y2), (x3, y3) = XY[list(c)] det = (x2 - x1) * (y3 - y1) - (x3 - x1) * (y2 - y1) assert det > 1e-12, "flipped or degenerate element; reduce amplitude" A = 0.5 * det b = np.array([y2 - y3, y3 - y1, y1 - y2]) / det cc = np.array([x3 - x2, x1 - x3, x2 - x1]) / det idx = list(c) K[np.ix_(idx, idx)] += A * (np.outer(b, b) + np.outer(cc, cc)) Mq[np.ix_(idx, idx)] += A * ME for jl, vd in enumerate(idx): Bx[idx, vd] += -b[jl] * A / 3.0 By[idx, vd] += -cc[jl] * A / 3.0 interior = np.where(~bnd)[0] Ki = K[np.ix_(interior, interior)] Kv = np.block([[Ki, np.zeros_like(Ki)], [np.zeros_like(Ki), Ki]]) B = np.hstack([Bx[:, interior], By[:, interior]]) return Kv, B, Mq, XY def kernel_dim(B): s = np.linalg.svd(B.T, compute_uv=False) return B.shape[0] - int(np.sum(s > 1e-10 * s[0])) def spectrum(Kv, B, Mq): S = B @ np.linalg.solve(Kv, B.T) L = np.linalg.cholesky(Mq) Li = np.linalg.inv(L) Sh = Li @ S @ Li.T lam = np.linalg.eigvalsh(0.5 * (Sh + Sh.T)) nz = int(np.sum(np.abs(lam) < 1e-10 * lam[-1])) beta_tilde = np.sqrt(max(lam[nz], 0.0)) return lam, nz, beta_tilde def patch_check(N=8): """Matrix-free verification (independent code path, right diagonal, structured): patch sums of grad(q_h) vanish at interior vertices.""" Kv, B, Mq, XY = build(N, "right") _, s, Vt = np.linalg.svd(B.T, full_matrices=True) r = int(np.sum(s > 1e-10 * s[0])) q = Vt.T[:, r + 2].reshape(N + 1, N + 1) # some kernel vector h = 1.0 / N worst = 0.0 for i in range(2, N - 1): for j in range(2, N - 1): g = np.array([ ((q[j-1, i] - q[j-1, i-1]) / h, (q[j, i] - q[j-1, i]) / h), ((q[j, i] - q[j, i-1]) / h, (q[j, i-1] - q[j-1, i-1]) / h), ((q[j, i+1] - q[j, i]) / h, (q[j, i] - q[j-1, i]) / h), ((q[j, i] - q[j, i-1]) / h, (q[j+1, i] - q[j, i]) / h), ((q[j, i+1] - q[j, i]) / h, (q[j+1, i+1] - q[j, i+1]) / h), ((q[j+1, i+1] - q[j+1, i]) / h, (q[j+1, i] - q[j, i]) / h), ]) worst = max(worst, np.abs(g.sum(axis=0)).max()) return worst def main(): # F: counting violation at N = 4 Kv, B, Mq, _ = build(4) print(f"N=4 counting check: dim ker B^T = {kernel_dim(B)} (>= 7 forced)\n") # A + D: structured sweep, right diagonal Ns = [8, 12, 16, 24, 32, 40] print("structured (right diagonal):") print(" N h zero modes beta_tilde") bt = {} for N in Ns: Kv, B, Mq, _ = build(N, "right") lam, nz, b = spectrum(Kv, B, Mq) assert nz == 8, (N, nz) # the constant pressure is in ker B^T exactly (all velocity dofs # are interior, so each column of B sums to a vanishing boundary # integral): assert np.linalg.norm(np.ones(Mq.shape[0]) @ B) < 1e-12 bt[N] = b print(f" {N:3d} {1.0/N: .5f} {nz} {b: .6e}") logh = np.log([1.0 / N for N in Ns]) logb = np.log([bt[N] for N in Ns]) print(f" fitted decay rate, all points : {np.polyfit(logh, logb, 1)[0]:.3f}") print(f" fitted decay rate, last three : " f"{np.polyfit(logh[-3:], logb[-3:], 1)[0]:.3f}") assert all(bt[Ns[i]] > bt[Ns[i + 1]] for i in range(len(Ns) - 1)) assert bt[32] / bt[8] < 0.35 # A: diagonal agreement (mirror symmetry) for N in [8, 16]: Kv, B, Mq, _ = build(N, "left") lam, nz, bL = spectrum(Kv, B, Mq) assert nz == 8 and abs(bL - bt[N]) < 1e-9, (N, nz, bL, bt[N]) print(" left diagonal: identical (8 zeros; beta_tilde agrees to 1e-9)\n") # B: perturbed meshes -> 5 exact modes print("perturbed meshes (amplitude 0.2 h):") for N, amp_bnd, seeds in [(8, 0.0, [1, 2, 3]), (16, 0.0, [1]), (8, 0.2, [1, 2, 3])]: for sd in seeds: Kv, B, Mq, _ = build(N, "right", amp_int=0.2, amp_bnd=amp_bnd, seed=sd) k = kernel_dim(B) assert k == 5, (N, amp_bnd, sd, k) where = "interior only" if amp_bnd == 0 else "all vertices " print(f" N={N:2d}, {where}, seeds {seeds}: dim ker B^T = 5") print() # C: alternating pattern print("alternating-diagonal pattern:") for N in [8, 12]: _, B, _, _ = build(N, "alt") ks = kernel_dim(B) _, Bp, _, _ = build(N, "alt", amp_int=0.2, seed=0) kp = kernel_dim(Bp) assert ks == 8 and kp == 7, (N, ks, kp) print(f" N={N:2d}: structured {ks}, perturbed {kp}") print() # E: independent matrix-free check w = patch_check(8) print(f"matrix-free patch-sum check on a kernel vector: max = {w:.2e}") assert w < 1e-12 print("\nall assertions passed") if __name__ == "__main__": main()