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