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`. `solvers/richardson.m`'s `for k = 1:niter` loop *is*
68this: `Un^(0)` is the divide computed just before the loop, each pass
69computes `Un^(k+1)` from `Un^(k)`, and `dlap` — evaluated once per iteration —
70is its own function, `lib/dlap.m`. A model's step calls the solver once per
71species (`Un = solve(Bu, dt * D1, ...)`, routed to the selected solvers/*.m
72file by a host-generated shim), which is where the split pays: a different
73solver for the same operator is a selector change — or a different call in
74the model — with `lib/dlap.m` untouched. The solver is written as a full re-evaluation
75rather than an accumulated correction `δ = Un^(k+1) - Un^(k)` on purpose:
76where `dlap` evaluates to zero exactly, every `Un^(k)` is bit-for-bit
77`Un^(0)`, with no cancellation to round differently. (In practice `dlap` is a
78real computation through chained fp32 transforms, so on the round sphere it
79lands near zero rather than at it — the tests bound how near.)
81## Convergence, and why it isn't GMRES
83Writing `M = I - dt*D*lap_s` and `A = M - dt*D*dlap`, each step is
84`Un^(k+1) = M^{-1}(B + dt*D*dlap(Un^(k)))` — a stationary iteration that
85converges to the exact solution of `A·Un = B` exactly when the spectral
86radius of `M^{-1}(dt*D*dlap)` is below 1: while the geometric correction
87stays small against what the round-sphere solve already inverts. Unlike
88GMRES, there is no residual check and no adaptive iteration count: `niter` is
89fixed before the run starts, so a shape/timestep/diffusivity combination
90outside the convergence radius fails silently — the state saturates or
91diverges over many steps — rather than being caught the way algos.tex's
92`solve_step` catches it (its `info != 0` return, logged when GMRES fails to
93reach `tol` within `maxiter`).
95That tradeoff is deliberate, not an oversight, and it comes from where the
96two projects run. algos.tex's GMRES needs, every iteration: a dot product
97across the whole spectral state (Arnoldi orthogonalization) and a residual
98norm to test against `tol` — both require reading a scalar back to the host
99mid-solve. This project's solver instead records one whole timestep as a
100single GPU command buffer, submitted once, with the entire `for k = 1:niter`
101loop unrolled at compile time into a fixed sequence of dispatches — there is
102no point in that sequence where the host makes a decision, and no path for a
103data-dependent stopping rule to plug in. (Recompiling — which changing
104`niter` triggers — is the only way this project can change how much work a
105step does; see the README's "`for` loops, unrolled".) Richardson iteration is
106the cheapest method that still fits that shape: the same preconditioner as
107algos.tex, one `dlap` evaluation per iteration, a fixed and
108recompile-on-change trip count, in exchange for linear rather than
109superlinear convergence.
111A Krylov method fits the shape too, as long as its scalars stay on the GPU —
112which is what `solvers/bicgstab.m` does: `dot` is a reduction dispatch, the
113alpha/omega/rho recurrences are 1-element kernels, and the iteration count is
114still fixed and unrolled. What it cannot have is exactly what GMRES's `tol`
115gives algos.tex: a stopping rule. It compensates in two ways — every ratio is
116algebraically guarded (`a*b/(b*b + eps)`) so a converged iteration goes
117stationary rather than dividing noise by noise, and the cost is fixed at two
118`dlap` evaluations plus three reductions per iteration whether or not it has
119already converged. In exchange it converges superlinearly, including on
120shape/timestep combinations outside the Richardson iteration's spectral
121radius — the tests pin Schnakenberg on the peanut at the app's default lmax
122as exactly such a case.
124algos.tex's own method is here too, in the same fixed-count form:
125`solvers/gmres.m` is right-preconditioned GMRES(niter) — one Arnoldi sweep,
126Givens rotations, back-substitution — minus the restart loop and minus
127`tol`/`maxiter`, since there is still no data-dependent stopping. Its basis
128and Hessenberg bookkeeping run through the indexed-access ops
129(`getslab`/`setslab`, `getat`/`setat`), which the planner compiles to
130static-offset buffer copies once the unrolled loop's variable makes every
131index a literal; the same guarded-ratio discipline covers the rotation and
132back-substitution divides. Where algos.tex's GMRES stops at `tol`, this one
133spends its fixed niter·(niter+3)/2 reductions and niter `dlap` evaluations
134and keeps whatever residual that bought.
136## One more difference worth flagging
138algos.tex maps the half-spectrum complex coefficients through a real vector
139space embedding `E` (Sec. 6.4) because GMRES needs one flat, real-linear
140operator to hand to a generic solver. This project never needs `E`/`E^{-1}`:
141its spectral state is *already* carried as a real "2 x nlm" array — row 0 the
142real part, row 1 the imaginary — rather than packed complex, so every step
143here, `dlap` included, is already ℝ-linear arithmetic on that layout with no
144embedding or un-embedding step at all.