{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "f77cc0e9",
   "metadata": {
    "tags": [
     "header"
    ]
   },
   "source": [
    "# Lecture 4, notebook L4A — stable Stokes pairs: rates, the pencil, and the $\\varepsilon$-shift\n",
    "\n",
    "Mixed finite element methods — a crash course. Companion to slides 22–23\n",
    "(Demonstrations 1a, 1b, 1c).\n",
    "\n",
    "**Conventions of this course.**\n",
    "1. Cells are run top to bottom; every cell is self-contained given the setup cell.\n",
    "2. The first Firedrake solve in a session includes JIT compilation; quoted timings are from second solves.\n",
    "3. If Firedrake is unavailable, the fallback figures in `figs/` show the rehearsal output.\n",
    "4. All triangular meshes use the right-diagonal pattern\n",
    "(`diagonal=\"right\"`), the convention of the reference implementation.\n",
    "Firedrake's default is `\"left\"`; matrix quantities (the pencil, the\n",
    "spectra) are identical across the two diagonal families, but error\n",
    "constants are not, so the convention matters for reproducing values.\n",
    "5. Values marked **[verified]** were produced by the pure NumPy reference implementation `verify_stokes_2d.py`; the Firedrake runs below should reproduce them to the digits shown, except where a note states otherwise.\n",
    "\n",
    "**Note on Demonstration 1c (2026-07-11).** The demonstration was\n",
    "redesigned after verification: on an enclosed-flow problem with $g = 0$,\n",
    "interior forcing cannot excite the pressure kernel (the Schur right-hand\n",
    "side lies in $\\operatorname{range}(\\mathsf B) = (\\ker \\mathsf B^T)^\\perp$\n",
    "by construction). The cell below shows the surprise, the failure that\n",
    "*is* present, and the Lecture-1 mechanism reproduced on demand.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "0d8f8e37",
   "metadata": {
    "tags": [
     "setup"
    ]
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Firedrake 2026.4.1 (PETSc 3.25.0) ready; figures directory: figs/\n"
     ]
    }
   ],
   "source": [
    "# Setup 1/2: Firedrake (FEM-on-Colab guard) and imports.\n",
    "try:\n",
    "    import firedrake\n",
    "except ImportError:\n",
    "    import subprocess, urllib.request\n",
    "    urllib.request.urlretrieve(\n",
    "        \"https://fem-on-colab.github.io/releases/firedrake-install-release-real.sh\",\n",
    "        \"/tmp/firedrake-install.sh\")\n",
    "    subprocess.run([\"bash\", \"/tmp/firedrake-install.sh\"], check=True)\n",
    "    import firedrake\n",
    "\n",
    "from firedrake import *\n",
    "import numpy as np\n",
    "import scipy.linalg as sla\n",
    "import scipy.sparse as sp\n",
    "import scipy.sparse.linalg as spsla\n",
    "import matplotlib.pyplot as plt\n",
    "import os\n",
    "os.makedirs(\"figs\", exist_ok=True)\n",
    "import logging\n",
    "# The manufactured data has polynomial degree up to twelve; against P1\n",
    "# arguments TSFC warns that its (correct) quadrature estimate exceeds\n",
    "# the argument degree tenfold.  The estimate is what exact integration\n",
    "# requires, so the warning is cosmetic and silenced here.\n",
    "logging.getLogger(\"tsfc\").setLevel(logging.ERROR)\n",
    "try:\n",
    "    _fd_version = firedrake.__version__\n",
    "except AttributeError:\n",
    "    try:\n",
    "        from importlib.metadata import version as _pkg_version\n",
    "        _fd_version = _pkg_version(\"firedrake\")\n",
    "    except Exception:\n",
    "        _fd_version = \"unknown\"\n",
    "from firedrake.petsc import PETSc as _PETSc\n",
    "_petsc = \".\".join(map(str, _PETSc.Sys.getVersion()))\n",
    "print(f\"Firedrake {_fd_version} (PETSc {_petsc}) ready; \"\n",
    "      \"figures directory: figs/\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "1c47e03a",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "definitions loaded: nu, structured_unit_square (diagonal='right'), manufactured, perturbed_unit_square, velocity_space, pressure_space, stokes_solve, stokes_errors, rate_report\n"
     ]
    }
   ],
   "source": [
    "# Setup 2/2: common definitions.\n",
    "nu = Constant(1.0)   # kept as a symbol throughout the lecture\n",
    "\n",
    "def structured_unit_square(N):\n",
    "    \"\"\"The reference-implementation triangulation: right diagonals.\n",
    "    Firedrake's default is diagonal='left'; matrix quantities (the\n",
    "    pencil) are diagonal-independent, error constants are not.\"\"\"\n",
    "    return UnitSquareMesh(N, N, diagonal=\"right\")\n",
    "\n",
    "def manufactured(mesh):\n",
    "    \"\"\"Decision D4: psi = x^2 (1-x)^2 y^2 (1-y)^2, u = curl psi,\n",
    "    p = x^3 y^3 - 1/16 (zero mean), nu = 1.\"\"\"\n",
    "    x, y = SpatialCoordinate(mesh)\n",
    "    g  = lambda t: t**2 * (1 - t)**2\n",
    "    gp = lambda t: 2*t*(1 - t)*(1 - 2*t)\n",
    "    u_ex = as_vector([g(x) * gp(y), -gp(x) * g(y)])\n",
    "    p_ex = x**3 * y**3 - Constant(1.0 / 16.0)\n",
    "    f = -nu * div(grad(u_ex)) + grad(p_ex)\n",
    "    return u_ex, p_ex, f\n",
    "\n",
    "def perturbed_unit_square(N, amp=0.2, seed=20260711):\n",
    "    \"\"\"The reference-implementation perturbed mesh, reproduced exactly.\n",
    "\n",
    "    The random stream, the amplitude, and the seed match\n",
    "    verify_stokes_2d.tri_mesh; the shifts are generated over the interior\n",
    "    vertices in the reference (lexicographic) ordering and then assigned\n",
    "    to Firedrake's vertices through the grid index, because the two\n",
    "    implementations enumerate vertices differently.  Same realization,\n",
    "    hence digit-level reproduction of the perturbed-family records.\"\"\"\n",
    "    mesh = structured_unit_square(N)\n",
    "    xs = np.linspace(0.0, 1.0, N + 1)\n",
    "    pts = np.stack(np.meshgrid(xs, xs, indexing=\"ij\"), axis=-1).reshape(-1, 2)\n",
    "    rng = np.random.default_rng(seed)\n",
    "    interior = (pts[:, 0] > 1e-12) & (pts[:, 0] < 1 - 1e-12) \\\n",
    "             & (pts[:, 1] > 1e-12) & (pts[:, 1] < 1 - 1e-12)\n",
    "    pts[interior] += (rng.random((interior.sum(), 2)) - 0.5) * 2 * amp / N\n",
    "    X = mesh.coordinates.dat.data\n",
    "    idx = np.rint(X[:, 0] * N).astype(int) * (N + 1) \\\n",
    "        + np.rint(X[:, 1] * N).astype(int)\n",
    "    X[:] = pts[idx]\n",
    "    return mesh\n",
    "\n",
    "def velocity_space(mesh, pair):\n",
    "    cell = mesh.ufl_cell()\n",
    "    if pair == \"TH\":\n",
    "        return VectorFunctionSpace(mesh, \"CG\", 2)\n",
    "    if pair == \"P2P0\":            # velocity as TH; the pressure differs\n",
    "        return VectorFunctionSpace(mesh, \"CG\", 2)\n",
    "    if pair == \"MINI\":\n",
    "        P1 = FiniteElement(\"CG\", cell, 1)\n",
    "        B  = FiniteElement(\"B\",  cell, 3)          # one bubble per element\n",
    "        return VectorFunctionSpace(mesh, P1 + B)\n",
    "    if pair == \"P1\":\n",
    "        return VectorFunctionSpace(mesh, \"CG\", 1)\n",
    "    raise ValueError(pair)\n",
    "\n",
    "def pressure_space(mesh, pair):\n",
    "    if pair == \"P2P0\":\n",
    "        return FunctionSpace(mesh, \"DG\", 0)   # piecewise constants\n",
    "    return FunctionSpace(mesh, \"CG\", 1)\n",
    "\n",
    "def stokes_solve(mesh, pair):\n",
    "    \"\"\"Direct (LU) solve of the manufactured Stokes problem; the constant\n",
    "    pressure is removed through the nullspace, not by pinning a dof.\"\"\"\n",
    "    W = velocity_space(mesh, pair) * pressure_space(mesh, pair)\n",
    "    u_ex, p_ex, f = manufactured(mesh)\n",
    "    w = Function(W)\n",
    "    u, p = TrialFunctions(W); v, q = TestFunctions(W)\n",
    "    a = nu*inner(grad(u), grad(v))*dx - p*div(v)*dx - q*div(u)*dx\n",
    "    L = inner(f, v)*dx\n",
    "    bcs = DirichletBC(W.sub(0), Constant((0.0, 0.0)), \"on_boundary\")\n",
    "    ns = MixedVectorSpaceBasis(\n",
    "        W, [W.sub(0), VectorSpaceBasis(constant=True, comm=W.comm)])\n",
    "    solve(a == L, w, bcs=bcs, nullspace=ns,\n",
    "          solver_parameters={\"mat_type\": \"aij\", \"ksp_type\": \"preonly\",\n",
    "                             \"pc_type\": \"lu\",\n",
    "                             \"pc_factor_mat_solver_type\": \"mumps\"})\n",
    "    return w\n",
    "\n",
    "def stokes_errors(w, mesh):\n",
    "    u_ex, p_ex, _ = manufactured(mesh)\n",
    "    uh, ph = w.subfunctions\n",
    "    pbar = assemble(ph * dx)                       # |Omega| = 1\n",
    "    e_grad = sqrt(assemble(inner(grad(uh - u_ex), grad(uh - u_ex)) * dx))\n",
    "    e_p    = sqrt(assemble((ph - pbar - p_ex)**2 * dx))\n",
    "    e_u    = sqrt(assemble(inner(uh - u_ex, uh - u_ex) * dx))\n",
    "    return e_grad, e_p, e_u\n",
    "\n",
    "def rate_report(errs, Ns, label):\n",
    "    errs = np.array(errs)\n",
    "    print(f\"{label}:  N    |grad(u-uh)|   |p-ph|      |u-uh|\")\n",
    "    for k, N in enumerate(Ns):\n",
    "        tail = \"\"\n",
    "        if k > 0:\n",
    "            r = np.log2(errs[k-1] / errs[k])\n",
    "            tail = f\"   rates: {r[0]:.2f}, {r[1]:.2f}, {r[2]:.2f}\"\n",
    "        print(f\"      {N:3d}  {errs[k,0]:.4e}  {errs[k,1]:.4e}\"\n",
    "              f\"  {errs[k,2]:.4e}{tail}\")\n",
    "\n",
    "print(\"definitions loaded: nu, structured_unit_square (diagonal='right'), \"\n",
    "      \"manufactured, perturbed_unit_square, velocity_space, \"\n",
    "      \"pressure_space, stokes_solve, stokes_errors, rate_report\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "32a44d30",
   "metadata": {
    "tags": [
     "demo1a"
    ]
   },
   "source": [
    "## Demonstration 1a (slide 22) — convergence rates\n",
    "\n",
    "The three stable pairs on the manufactured solution, structured meshes,\n",
    "$N = 8$ to $64$. Every rate in the table is a statement from the lecture:\n",
    "the $P_2$–$P_0$ row is Theorem 35, the MINI row is Corollary 39\n",
    "together with the duality remark, the Taylor–Hood row is Theorem 42.\n",
    "\n",
    "**Expected output [verified]** (last-step rates, structured meshes):\n",
    "\n",
    "| pair | $\\|\\nabla(u-u_h)\\|_0$ | $\\|p-p_h\\|_0$ | $\\|u-u_h\\|_0$ |\n",
    "|---|---|---|---|\n",
    "| $P_2$–$P_0$ | 0.97 | 1.01 | 1.94 |\n",
    "| MINI | 1.01 | 1.53 | 2.01 |\n",
    "| Taylor–Hood | 2.00 | 2.00 | 3.00 |\n",
    "\n",
    "The MINI pressure rate exceeds the guaranteed first order on these\n",
    "structured meshes (measured 1.53; on perturbed meshes, 1.27 — see the\n",
    "self-study cell). Corollary 39 is a guarantee from below; the observed\n",
    "extra half order is reported as measured. The $P_2$–$P_0$ rates are\n",
    "first order in energy and pressure, limited by the piecewise constant\n",
    "pressure (Lecture 2, Corollary 8).\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "33071c67",
   "metadata": {
    "tags": [
     "demo1a"
    ]
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "P2P0:  N    |grad(u-uh)|   |p-ph|      |u-uh|\n",
      "        8  2.1638e-02  2.7134e-02  8.0605e-04\n",
      "       16  1.1813e-02  1.3339e-02  2.3462e-04   rates: 0.87, 1.02, 1.78\n",
      "       32  6.1760e-03  6.5534e-03  6.3411e-05   rates: 0.94, 1.03, 1.89\n",
      "       64  3.1557e-03  3.2434e-03  1.6477e-05   rates: 0.97, 1.01, 1.94\n",
      "\n",
      "MINI:  N    |grad(u-uh)|   |p-ph|      |u-uh|\n",
      "        8  1.9044e-02  1.1653e-02  8.8791e-04\n",
      "       16  9.4877e-03  3.9037e-03  2.2333e-04   rates: 1.01, 1.58, 1.99\n",
      "       32  4.7123e-03  1.3129e-03  5.5280e-05   rates: 1.01, 1.57, 2.01\n",
      "       64  2.3465e-03  4.5449e-04  1.3719e-05   rates: 1.01, 1.53, 2.01\n",
      "\n",
      "TH:  N    |grad(u-uh)|   |p-ph|      |u-uh|\n",
      "        8  2.6145e-03  2.7160e-03  4.3887e-05\n",
      "       16  6.5746e-04  6.8341e-04  5.3449e-06   rates: 1.99, 1.99, 3.04\n",
      "       32  1.6461e-04  1.7117e-04  6.6389e-07   rates: 2.00, 2.00, 3.01\n",
      "       64  4.1170e-05  4.2813e-05  8.2876e-08   rates: 2.00, 2.00, 3.00\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# Demo 1a: rates on structured meshes.\n",
    "Ns = [8, 16, 32, 64]\n",
    "for pair in (\"P2P0\", \"MINI\", \"TH\"):\n",
    "    errs = []\n",
    "    for N in Ns:\n",
    "        mesh = structured_unit_square(N)\n",
    "        errs.append(stokes_errors(stokes_solve(mesh, pair), mesh))\n",
    "    rate_report(errs, Ns, pair)\n",
    "    print()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ba88de99",
   "metadata": {
    "tags": [
     "demo1b"
    ]
   },
   "source": [
    "## Demonstration 1b (slide 23) — the pencil, fourth appearance\n",
    "\n",
    "The generalized eigenvalue problem\n",
    "$\\mathsf B \\mathsf K_V^{-1} \\mathsf B^T \\mathbf q = \\lambda \\mathsf M_p \\mathbf q$\n",
    "on zero-mean pressures (Lecture 2, Exercise 1(a)), measured for the\n",
    "Stokes pairs; $\\beta_h = \\sqrt{\\lambda_{\\min}}$.\n",
    "\n",
    "**Expected output [verified]** ($N = 8, 16, 24$; structured meshes):\n",
    "MINI $\\beta_h = 0.3143,\\ 0.3136,\\ 0.3134$; Taylor–Hood\n",
    "$\\beta_h = 0.3662,\\ 0.3656,\\ 0.3654$ — flat under refinement, which is\n",
    "the assertion. For $P_1$–$P_1$ at $N = 8$: $\\dim\\ker\\mathsf B^T = 8$\n",
    "including the constant, that is, $7$ modes beyond it (the live assertion\n",
    "of slide 23). The largest pencil eigenvalue stays below $1$, which is\n",
    "Lemma 43 in numerical form (Theorem 45; used again in `L4B`).\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "d34f0632",
   "metadata": {
    "tags": [
     "demo1b"
    ]
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "MINI: beta_h = 0.3143, 0.3136, 0.3134   (N = 8, 16, 24)\n",
      "TH  : beta_h = 0.3662, 0.3656, 0.3654   (N = 8, 16, 24)\n",
      "P1-P1, N = 8: dim ker B^T = 8 (including the constant), that is, 7 beyond it\n"
     ]
    }
   ],
   "source": [
    "# Demo 1b: the pencil.  Matrices are assembled on separate (non-mixed)\n",
    "# spaces, exported to SciPy, and restricted to interior velocity dofs,\n",
    "# exactly as in the reference implementation.\n",
    "\n",
    "def petsc_to_csr(A):\n",
    "    ai, aj, av = A.petscmat.getValuesCSR()\n",
    "    return sp.csr_matrix((av, aj, ai), shape=A.petscmat.getSize())\n",
    "\n",
    "def pencil(mesh, pair):\n",
    "    V = velocity_space(mesh, pair)\n",
    "    Q = FunctionSpace(mesh, \"CG\", 1)\n",
    "    u, v = TrialFunction(V), TestFunction(V)\n",
    "    p, q = TrialFunction(Q), TestFunction(Q)\n",
    "    KV = petsc_to_csr(assemble(inner(grad(u), grad(v))*dx, mat_type=\"aij\"))\n",
    "    B  = petsc_to_csr(assemble(-div(u)*q*dx, mat_type=\"aij\"))   # rows: Q\n",
    "    Mp = petsc_to_csr(assemble(p*q*dx, mat_type=\"aij\"))\n",
    "    bd = DirichletBC(V, 0, \"on_boundary\").nodes\n",
    "    bs = V.value_size\n",
    "    mask = np.ones(V.dim(), bool)\n",
    "    for c in range(bs):\n",
    "        mask[bs*np.asarray(bd) + c] = False\n",
    "    iv = np.nonzero(mask)[0]\n",
    "    KVi = KV[iv][:, iv].tocsc(); Bi = B[:, iv].tocsr()\n",
    "    S = Bi @ spsla.splu(KVi).solve(Bi.T.toarray())\n",
    "    S = 0.5*(S + S.T)\n",
    "    m = np.asarray(Mp @ np.ones(Mp.shape[0])).ravel()\n",
    "    Z = sla.null_space(m[None, :])                 # zero-mean pressures\n",
    "    lam = sla.eigh(Z.T @ S @ Z, Z.T @ (Mp @ Z), eigvals_only=True)\n",
    "    sv = sla.svdvals(Bi.toarray())\n",
    "    kdim = Bi.shape[0] - int((sv > 1e-9*sv.max()).sum())\n",
    "    return np.sort(lam), kdim\n",
    "\n",
    "for pair in (\"MINI\", \"TH\"):\n",
    "    vals = []\n",
    "    for N in (8, 16, 24):\n",
    "        lam, kdim = pencil(structured_unit_square(N), pair)\n",
    "        assert kdim == 1                            # the constant only\n",
    "        vals.append(np.sqrt(lam[0]))\n",
    "    print(f\"{pair:4s}: beta_h =\", \", \".join(f\"{v:.4f}\" for v in vals),\n",
    "          \"  (N = 8, 16, 24)\")\n",
    "\n",
    "lam, kdim = pencil(structured_unit_square(8), \"P1\")\n",
    "print(f\"P1-P1, N = 8: dim ker B^T = {kdim} (including the constant), \"\n",
    "      f\"that is, {kdim - 1} beyond it\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3ca8e0b4",
   "metadata": {
    "tags": [
     "demo1c"
    ]
   },
   "source": [
    "## Demonstration 1c (slide 23, revised) — the $\\varepsilon$-shift, reread\n",
    "\n",
    "Three measurements.\n",
    "\n",
    "**(i) The surprise.** The $\\varepsilon$-shifted $P_1$–$P_1$ solve of the\n",
    "manufactured problem produces a *clean-looking* pressure:\n",
    "$\\max|p_h| \\approx 1.006$ **[verified]**, independent of $\\varepsilon$.\n",
    "Reason: with $g = 0$ the Schur right-hand side\n",
    "$\\mathsf B \\mathsf K^{-1}\\mathbf f$ lies in $\\operatorname{range}(\\mathsf B)\n",
    "= (\\ker \\mathsf B^T)^\\perp$ *by construction* — interior forcing can never\n",
    "excite the kernel. The Lecture-1 cavity amplitude $763$ entered through\n",
    "the inhomogeneous lid lifting in the constraint, not through $\\mathbf f$.\n",
    "\n",
    "**(ii) The failure that is present.** The pair is still singular, and the\n",
    "essential inf-sup constant still decays. On the manufactured problem this\n",
    "appears as *pressure stagnation*: **[verified]**, structured meshes,\n",
    "$\\|p - p_h\\|_0 \\approx 1.06,\\ 0.88,\\ 0.82,\\ 0.80 \\times 10^{-1}$ for\n",
    "$N = 8, 16, 32, 64$ (rates $0.26, 0.11, 0.05$), while the velocity\n",
    "converges at the healthy rates $1.01$ (energy) and $2.01$ ($L^2$).\n",
    "The pressure does not converge; nothing in the picture warns of it.\n",
    "\n",
    "**(iii) The mechanism, on demand.** Injecting a kernel component of size\n",
    "$|g_k| = 10^{-3}$ into the constraint datum reproduces the Lecture-1\n",
    "blow-up: amplitude of order $10^8$ for $\\varepsilon = 10^{-8}$\n",
    "(reference implementation: $3.288\\cdot 10^8$; the exact value depends on\n",
    "which unit kernel mode the SVD returns, so only the order and the scaling\n",
    "$\\propto |g_k|/\\varepsilon$ are asserted here — verify the scaling by\n",
    "changing $\\delta$ or $\\varepsilon$ by a factor of ten).\n",
    "\n",
    "The figure under (i) shows the point precisely: the pressure is\n",
    "*bounded* (the amplitude matches $\\max|p^*|$), which is the surprise,\n",
    "but it is visibly polluted at the mesh scale — that oscillatory\n",
    "component is exactly the non-convergent part quantified by the\n",
    "stagnating pressure error in (ii).\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e2169deb",
   "metadata": {
    "tags": [
     "demo1c"
    ]
   },
   "outputs": [],
   "source": [
    "# Demo 1c: (i) shifted solve, (ii) stagnation, (iii) kernel injection.\n",
    "\n",
    "def p1p1_shift(mesh, eps=1e-8, g_dual=None):\n",
    "    \"\"\"eps-shifted P1-P1 solve; g_dual, if given, is added to the\n",
    "    assembled pressure right-hand side (a raw dual vector, matching the\n",
    "    reference implementation).\"\"\"\n",
    "    W = velocity_space(mesh, \"P1\") * FunctionSpace(mesh, \"CG\", 1)\n",
    "    u_ex, p_ex, f = manufactured(mesh)\n",
    "    w = Function(W)\n",
    "    u, p = TrialFunctions(W); v, q = TestFunctions(W)\n",
    "    a = nu*inner(grad(u), grad(v))*dx - p*div(v)*dx - q*div(u)*dx \\\n",
    "        - Constant(eps)*p*q*dx\n",
    "    L = inner(f, v)*dx\n",
    "    bcs = DirichletBC(W.sub(0), Constant((0.0, 0.0)), \"on_boundary\")\n",
    "    A = assemble(a, bcs=bcs, mat_type=\"aij\")\n",
    "    b = assemble(L, bcs=bcs)\n",
    "    if g_dual is not None:\n",
    "        b.subfunctions[1].dat.data[:] += g_dual\n",
    "    solve(A, w, b, solver_parameters={\"ksp_type\": \"preonly\",\n",
    "                                      \"pc_type\": \"lu\"})\n",
    "    ph = w.subfunctions[1]\n",
    "    pbar = assemble(ph * dx)\n",
    "    return w, float(np.abs(ph.dat.data - pbar).max())\n",
    "\n",
    "# (i) the surprise\n",
    "mesh = structured_unit_square(32)\n",
    "w, amp = p1p1_shift(mesh)\n",
    "print(f\"(i)  eps-shift on the manufactured problem: max|p_h| = {amp:.3e}\")\n",
    "fig, ax = plt.subplots(figsize=(4.0, 3.2))\n",
    "from firedrake.pyplot import tripcolor\n",
    "c = tripcolor(w.subfunctions[1], axes=ax, cmap=\"RdBu_r\")\n",
    "fig.colorbar(c, ax=ax); ax.set_aspect(\"equal\")\n",
    "ax.set_title(\"P1-P1, eps-shift: bounded amplitude, mesh-scale pollution\")\n",
    "plt.show()\n",
    "\n",
    "# (ii) the stagnation\n",
    "errs, Ns = [], [8, 16, 32, 64]\n",
    "for N in Ns:\n",
    "    mesh = structured_unit_square(N)\n",
    "    w, _ = p1p1_shift(mesh)\n",
    "    errs.append(stokes_errors(w, mesh))\n",
    "rate_report(errs, Ns, \"(ii) P1-P1\")\n",
    "\n",
    "# (iii) the mechanism: one unit kernel mode into the constraint datum\n",
    "mesh = structured_unit_square(32)\n",
    "V = velocity_space(mesh, \"P1\"); Q = FunctionSpace(mesh, \"CG\", 1)\n",
    "u, v = TrialFunction(V), TestFunction(V); q = TestFunction(Q)\n",
    "B = petsc_to_csr(assemble(-div(u)*q*dx, mat_type=\"aij\"))\n",
    "bd = DirichletBC(V, 0, \"on_boundary\").nodes\n",
    "mask = np.ones(V.dim(), bool)\n",
    "for c in range(V.value_size):\n",
    "    mask[V.value_size*np.asarray(bd) + c] = False\n",
    "Bi = B[:, np.nonzero(mask)[0]]\n",
    "U_, sv, _ = sla.svd(Bi.toarray(), full_matrices=True)\n",
    "qk = U_[:, int((sv > 1e-9*sv.max()).sum())]        # one unit kernel mode\n",
    "w, amp = p1p1_shift(mesh, g_dual=1e-3 * qk)\n",
    "print(f\"(iii) |g_k| = 1e-3, eps = 1e-8: max|p_h| = {amp:.3e}\"\n",
    "      \"   (order 1e8; scaling ~ |g_k|/eps)\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "71afea01",
   "metadata": {
    "tags": [
     "selfstudy"
    ]
   },
   "source": [
    "## Self-study\n",
    "\n",
    "**(a) Rates on perturbed meshes.** Rerun Demonstration 1a with\n",
    "`perturbed_unit_square(N)` in place of `structured_unit_square(N)`.\n",
    "**Expected [verified]** (last-step rates): $P_2$–$P_0$\n",
    "$0.97 / 1.01 / 1.95$, MINI $1.01 / 1.27 / 2.01$, Taylor–Hood\n",
    "$1.98 / 1.99 / 2.97$. The MINI pressure keeps part of its\n",
    "extra half order; the guarantee of Corollary 39 is first order.\n",
    "\n",
    "**(b) $Q_1$–$P_0$, measurement only.** On quadrilateral meshes\n",
    "(`UnitSquareMesh(N, N, quadrilateral=True)`, spaces `(\"CG\", 1)` and\n",
    "`(\"DQ\", 0)`), the pencil has *two* exact zero modes — the constant and\n",
    "the checkerboard (Lecture 1). Remove both and track the smallest\n",
    "remaining eigenvalue under refinement. No expected value is supplied;\n",
    "report what you measure. (Exercise sheet 4, problem 4, develops this.)\n",
    "\n",
    "**(c) Gradient forcing.** For $\\boldsymbol f = \\nabla\\phi$ with\n",
    "$\\phi = x^3 y^3 - \\tfrac{1}{16}$ the exact solution is\n",
    "$(\\boldsymbol u, p) = (\\boldsymbol 0, \\phi)$; Proposition 44 bounds the\n",
    "discrete velocity by $\\nu^{-1}\\inf_{q_h}\\|\\phi - q_h\\|_0$. The cell\n",
    "below prints the measured pollution $\\|\\nabla \\boldsymbol u_h\\|_0$, the\n",
    "projection bound, and their ratio on perturbed meshes.\n",
    "**Expected [verified]** ($\\nu = 1$; [V6]): at $N = 16, 32, 64$ the\n",
    "pollution is $4.95 / 1.31 / 0.336 \\times 10^{-4}$ for MINI and\n",
    "$2.27 / 0.694 / 0.184 \\times 10^{-4}$ for Taylor–Hood, against the\n",
    "bound $7.43 / 1.92 / 0.485 \\times 10^{-4}$; all columns converge with\n",
    "order two, and the ratio stays below one (slide 26).\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "14749e85",
   "metadata": {
    "tags": [
     "selfstudy"
    ]
   },
   "outputs": [],
   "source": [
    "# Self-study (a): perturbed-mesh rates -- complete code, run as is.\n",
    "for pair in (\"P2P0\", \"MINI\", \"TH\"):\n",
    "    errs = []\n",
    "    for N in Ns:\n",
    "        mesh = perturbed_unit_square(N)\n",
    "        errs.append(stokes_errors(stokes_solve(mesh, pair), mesh))\n",
    "    rate_report(errs, Ns, pair + \" (perturbed)\")\n",
    "    print()\n",
    "\n",
    "# Self-study (b): Q1-P0 essential constant -- measurement only.\n",
    "# Sketch: build the pencil on quadrilateral meshes; sla.eigh returns the\n",
    "# spectrum; discard the two smallest eigenvalues (the exact modes) and\n",
    "# tabulate the third under refinement.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9b40b7c7",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Self-study (c): gradient forcing -- complete code, run as is.\n",
    "def gradient_forcing(mesh, pair):\n",
    "    W = velocity_space(mesh, pair) * pressure_space(mesh, pair)\n",
    "    x, y = SpatialCoordinate(mesh)\n",
    "    phi = x**3 * y**3 - Constant(1.0/16.0)\n",
    "    w = Function(W)\n",
    "    u, p = TrialFunctions(W); v, q = TestFunctions(W)\n",
    "    a = nu*inner(grad(u), grad(v))*dx - p*div(v)*dx - q*div(u)*dx\n",
    "    L = inner(grad(phi), v)*dx           # f = grad(phi): u = 0, p = phi\n",
    "    bcs = DirichletBC(W.sub(0), Constant((0.0, 0.0)), \"on_boundary\")\n",
    "    ns = MixedVectorSpaceBasis(\n",
    "        W, [W.sub(0), VectorSpaceBasis(constant=True, comm=W.comm)])\n",
    "    solve(a == L, w, bcs=bcs, nullspace=ns,\n",
    "          solver_parameters={\"mat_type\": \"aij\", \"ksp_type\": \"preonly\",\n",
    "                             \"pc_type\": \"lu\",\n",
    "                             \"pc_factor_mat_solver_type\": \"mumps\"})\n",
    "    uh = w.subfunctions[0]\n",
    "    pollution = sqrt(assemble(inner(grad(uh), grad(uh))*dx))\n",
    "    qh = project(phi, pressure_space(mesh, pair))   # L2 projection\n",
    "    bound = sqrt(assemble((phi - qh)**2*dx)) / float(nu)\n",
    "    return pollution, bound\n",
    "\n",
    "print(\"gradient forcing, perturbed meshes, nu = 1   [verified: V6]\")\n",
    "print(\"   N   pair    |grad u_h|     bound          ratio\")\n",
    "for N in (16, 32, 64):\n",
    "    mesh = perturbed_unit_square(N)\n",
    "    for pair in (\"MINI\", \"TH\"):\n",
    "        pol, bnd = gradient_forcing(mesh, pair)\n",
    "        print(f\"  {N:3d}   {pair:4s}   {pol:.3e}     {bnd:.3e}\"\n",
    "              f\"     {pol/bnd:.3f}\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "85acc164-051f-43bf-8c58-d7b8cf45e84c",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Firedrake",
   "language": "python",
   "name": "firedrake"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.12.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}