1% Solve the implicit diffusion system
2%
3% (I - dtD*lap_g) X = B
4%
5% by preconditioned Richardson iteration, with the round-sphere operator as
6% the preconditioner. Splitting lap_g = lap_s + dlap and moving the geometric
7% part to the right-hand side gives the fixed point
8%
9% X = (B + dtD*dlap(X)) ./ (1 + dtD*lam)
10%
11% iterated from the round-sphere answer (lam holds +l(l+1), so the divide is
12% the exact inverse of I - dtD*lap_s). It converges while dtD*dlap stays
13% small against what the divide already inverts; niter is fixed at compile
14% time — the loop is unrolled into the op sequence, so there is no residual
15% check and no adaptive stopping. Full derivation and the map onto
16% algos.tex's GMRES formulation: docs/richardson-iteration.md.
17%
18% Written as a full re-evaluation rather than an accumulated correction on
19% purpose: on the round sphere dlap is identically zero, so every iterate is
20% bit for bit the first divide, with no cancellation to round differently.
22function X = richardson(B, dtD, lam, filt, Vtx, Vty, Vtz, Vpx, Vpy, Vpz, niter)
23 X = B ./ (1 + dtD * lam);
24 for k = 1:niter
25 dL = dlap(X, filt, Vtx, Vty, Vtz, Vpx, Vpy, Vpz, lam);
26 X = (B + dtD * dL) ./ (1 + dtD * lam);
27 end
28end