1# Notes on this project's Richardson iteration, for readers of `algos.tex`
3## Why this exists
5`evolving_surface/notes/algos.tex` (Sec. 5, "Implicit timestepping and the
6linear solve") specifies the surface diffusion step as backward Euler,
7`(I - Δt Δ_Γ) u^{n+1} = u^n`, solved by preconditioned GMRES: the
8round-sphere Laplacian `M = I - Δt Δ_S` (diagonal, invertible by eigenvalue)
9preconditions the full operator, the real-embedding map `E` puts the
10half-spectrum complex coefficients into a real vector space, and restarted
11GMRES iterates to a residual tolerance.
13turing-surface (this project) solves the *same* split operator with a
14different numerical method: a preconditioned Richardson (fixed-point)
15iteration, reusing exactly algos.tex's preconditioner `M^{-1}` but with no
16Krylov subspace, no orthogonalization, and no adaptive stopping. This note
17gives the map between the two, in this project's variable names, and why the
18switch.
20## Notation map
22| algos.tex | this project | meaning |
23|---|---|---|
24| `u^n`, `u^{n+1}` | `U`/`V` (in), `Un`/`Vn` (out) | spectral state, one array per species |
25| `Δ_Γ` | `lap_g` | the surface's Laplace-Beltrami operator |
26| `Δ_S` | `lap_s` | the round sphere's operator, eigenvalue `-l(l+1)` |
27| — | `dlap` | `lap_g - lap_s`. algos.tex has no name for this because it never splits the operator this way — its GMRES matvec (`surface_screened_laplacian`) applies the *whole* `Δ_Γ` every iteration. |
28| `M = I - Δt Δ_S` | `(1 + dt*D*lam)` | the same preconditioner. `lam` holds `+l(l+1)`, not `-l(l+1)`, so it enters as a *sum* — the sign flip is already folded into `lam`. |
29| `M^{-1}v` (eq. `preconditioner_inverse`) | `v ./ (1 + dt*D*lam)` | the identical elementwise divide |
30| a GMRES iterate | `Un^(k)`, `k = 0..niter` | *not* a Krylov iterate — a fresh, full re-solve of the fixed point below, evaluated at the previous iterate |
32## The fixed point this project actually iterates
34Same split as algos.tex, `lap_g = lap_s + dlap`, substituted into backward
35Euler and rearranged so every occurrence of the unknown is `Un`. Starting
36from `(I - dt*D*lap_g) Un = B` and substituting the split:
38```
39(I - dt*D*(lap_s + dlap)) Un = B
40```
42Expanding, and moving the `dlap` term to the right so only the exactly
43invertible round-sphere part remains on the left:
45```
46Un - dt*D*lap_s(Un) = B + dt*D*dlap(Un)
47```
49`lap_s` is diagonal with eigenvalue `-l(l+1)`, and `lam` holds `+l(l+1)`, so
50`lap_s(Un) = -lam .* Un` — the left side becomes `Un .* (1 + dt*D*lam)`, and
51dividing through gives:
53```
54Un = (B + dt*D*dlap(Un)) ./ (1 + dt*D*lam)
55```
57`B` is the explicit-reaction right-hand side — this project's models are
58IMEX (explicit reaction, implicit diffusion), where algos.tex's worked
59example is the bare heat equation, so `B` here is `u^n` plus a reaction term.
60Richardson iteration on this fixed point:
62```
63Un^(0) = B ./ (1 + dt*D*lam) [dlap = 0]
64Un^(k+1) = (B + dt*D*dlap(Un^(k))) ./ (1 + dt*D*lam)
65```
67for `k = 0 .. niter-1`. `models/schnakenberg.m`'s `for k = 1:niter` loop *is*
68this: `Un^(0)` is the divide computed just before the loop, and each pass
69computes `Un^(k+1)` from `Un^(k)`. It is written as a full re-evaluation
70rather than an accumulated correction `δ = Un^(k+1) - Un^(k)` on purpose: at
71`dlap ≡ 0` (the round sphere), every `Un^(k)` is then bit-for-bit `Un^(0)`,
72with no cancellation to round differently — a stronger, and cheaper to
73check, statement than "close to the round-sphere answer."
75## Convergence, and why it isn't GMRES
77Writing `M = I - dt*D*lap_s` and `A = M - dt*D*dlap`, each step is
78`Un^(k+1) = M^{-1}(B + dt*D*dlap(Un^(k)))` — a stationary iteration that
79converges to the exact solution of `A·Un = B` exactly when the spectral
80radius of `M^{-1}(dt*D*dlap)` is below 1: while the geometric correction
81stays small against what the round-sphere solve already inverts. Unlike
82GMRES, there is no residual check and no adaptive iteration count: `niter` is
83fixed before the run starts, so a shape/timestep/diffusivity combination
84outside the convergence radius fails silently — the state saturates or
85diverges over many steps — rather than being caught the way algos.tex's
86`solve_step` catches it (its `info != 0` return, logged when GMRES fails to
87reach `tol` within `maxiter`).
89That tradeoff is deliberate, not an oversight, and it comes from where the
90two projects run. algos.tex's GMRES needs, every iteration: a dot product
91across the whole spectral state (Arnoldi orthogonalization) and a residual
92norm to test against `tol` — both require reading a scalar back to the host
93mid-solve. This project's solver instead records one whole timestep as a
94single GPU command buffer, submitted once, with the entire `for k = 1:niter`
95loop unrolled at compile time into a fixed sequence of dispatches — there is
96no point in that sequence where the host makes a decision, and no path for a
97data-dependent stopping rule to plug in. (Recompiling — which changing
98`niter` triggers — is the only way this project can change how much work a
99step does; see the README's "`for` loops, unrolled".) Richardson iteration is
100the cheapest method that still fits that shape: the same preconditioner as
101algos.tex, one `dlap` evaluation per iteration, a fixed and
102recompile-on-change trip count, in exchange for linear rather than
103superlinear convergence.
105## One more difference worth flagging
107algos.tex maps the half-spectrum complex coefficients through a real vector
108space embedding `E` (Sec. 6.4) because GMRES needs one flat, real-linear
109operator to hand to a generic solver. This project never needs `E`/`E^{-1}`:
110its spectral state is *already* carried as a real "2 x nlm" array — row 0 the
111real part, row 1 the imaginary — rather than packed complex, so every step
112here, `dlap` included, is already ℝ-linear arithmetic on that layout with no
113embedding or un-embedding step at all.