0af3386Reaction-diffusion on spherical-harmonic surfacesJeremy Magland 1% Schnakenberg reaction-diffusion on a closed surface.
2%
3% du/dt = D1*lap_g(u) + a - u + u^2*v
4% dv/dt = D2*lap_g(v) + b - u^2*v
5%
62e6cc4Simplify UI text and built-in .m script commentsJeremy Magland 6% Explicit reaction, implicit diffusion (IMEX Euler). The implicit solve
7% splits lap_g = lap_s + dlap: the round-sphere part lap_s is diagonal in
8% spherical-harmonic space (eigenvalues -lam), and the loop iterates the
9% geometric correction dlap from that exact solve. Grid fields are npts x 1;
8144287WIP: added code for Laplace-Beltrami operator evaluation on smooth genus-0 surfaceOwen Melia 10% spectral fields are real 2 x nlm. See docs/richardson-iteration.md.
591a4f5Reduce the Laplace-Beltrami matvec to 6 transforms per species per iterationDan Fortunato 11%
12% The correction evaluates lap_g in flux form -- 6 transforms per species
13% per iteration where the Cartesian-gradient form (Algorithm 4 of
14% evolving_surface/notes/algos.tex) needs 12. See
15% docs/reduced-transforms.md, and models/schnakenberg_alg4.m
16% for the original form kept as a live reference.
0af3386Reaction-diffusion on spherical-harmonic surfacesJeremy Magland 17
18function [U, V, u, v] = init(noise, a, b)
19 us = a + b;
20 vs = b / (us * us);
a4fee9cBatch independent transforms through one Legendre dispatchDan Fortunato 21 [U, V] = analys(us + noise, vs * ones(numel(noise), 1));
22 [u, v] = synth(U, V);
e4d6a3bPrecondition with the operator's symbol; project the correction onto the bandDan Fortunato 25function [Un, Vn, u, v] = step(U, V, lam, filt, gx, gy, gz, p1, p2, q2, r, jhat, a, b, D1, D2, dt, niter)
a4fee9cBatch independent transforms through one Legendre dispatchDan Fortunato 26 % Grouped transforms -- [a, b] = synth(x, y) -- are explicit batching:
27 % output k is the transform of input k, and the whole group runs as one
28 % batched Legendre dispatch, or as many as the device's lane width allows
29 % (src/mgpu/plan.ts, materializeTransforms). The grouping is a promise of
30 % independence, never of a lane width, so the same source runs anywhere.
31 [u, v] = synth(U, V);
62e6cc4Simplify UI text and built-in .m script commentsJeremy Magland 34 % Right-hand side of the implicit solve (I - dt*D*lap_g) Unew = B.
36 rv = b - uuv;
37 [Ru, Rv] = analys(ru, rv);
38 Bu = U + dt * Ru;
39 Bv = V + dt * Rv;
e4d6a3bPrecondition with the operator's symbol; project the correction onto the bandDan Fortunato 41 % Preconditioned solve, then iterate the geometric correction. jhat is
42 % the host's minimax scale over the operator's symbol eigenvalues mu(x)
43 % -- the inverse squared principal stretches of the embedding, direction
44 % included (src/geom/geometry.ts, Jhat): preconditioning with lam/jhat
45 % contracts every mode and direction at rate
46 % (muMax - muMin)/(muMax + muMin) < 1 on any surface, where the plain
47 % lam (jhat = 1) diverges wherever mu > 2 -- docs/reduced-transforms.md
48 % Sec 10. The answer never depends on jhat (the lamJ term added inside
49 % dLu is the term divided back out); only the convergence rate does. On
50 % the sphere mu = 1 and lamJ = lam.
51 lamJ = lam ./ jhat;
52 Un = Bu ./ (1 + (dt * D1) * lamJ);
53 Vn = Bv ./ (1 + (dt * D2) * lamJ);
55 for k = 1:niter
591a4f5Reduce the Laplace-Beltrami matvec to 6 transforms per species per iterationDan Fortunato 56 % dlap = lap_g - lap_s at the current iterate, in flux form
a4fee9cBatch independent transforms through one Legendre dispatchDan Fortunato 57 % (docs/reduced-transforms.md Sec 4). The sin-weighted derivatives
58 % sin(theta)*dtheta(u) and dphi(u) -- both smooth on the sphere,
59 % synthesized straight from the dthetac/dphic coefficient shuffles --
60 % are combined pointwise through the precomputed weights p1,p2,q2 into
61 % two fluxes P,Q, also smooth. Their coefficients are then pushed
62 % through the *same* shuffles again and summed before the one synthesis
63 % of the divergence, which r scales into lap_g(u). The only division by
64 % sin(theta) anywhere is folded into p1,p2,q2,r at precompute time.
e4d6a3bPrecondition with the operator's symbol; project the correction onto the bandDan Fortunato 65 % lamJ.*Un adds back the preconditioner's -lap_s(Un)/jhat, since lam
66 % holds +l(l+1). filt zeroes the top two degrees, where the derivative
67 % recurrences cannot exactly represent a derivative -- and the correction
68 % itself is projected onto the same band (algos.tex Algorithm 5 zeroes
69 % the same coefficients): without that, each iteration replaces a bit
70 % more of the top degrees' implicit diffusion with nothing (their fixed
71 % point is the undiffused Bu), and the two species un-diffuse at
72 % different rates -- a spurious Turing band at the band edge.
74 % The two species share each grouped call: the four gradient
75 % syntheses, the four flux analyses, the two divergence syntheses and
76 % the two final analyses each run as one batched dispatch.
8144287WIP: added code for Laplace-Beltrami operator evaluation on smooth genus-0 surfaceOwen Melia 77 Fu = Un .* filt;
79 vtu = dthetac(Fu);
80 vpu = dphic(Fu);
81 vtv = dthetac(Fv);
82 vpv = dphic(Fv);
83 [Ftu, Fpu, Ftv, Fpv] = synth(vtu, vpu, vtv, vpv);
591a4f5Reduce the Laplace-Beltrami matvec to 6 transforms per species per iterationDan Fortunato 84 Pu = p1 .* Ftu + p2 .* Fpu;
85 Qu = p2 .* Ftu + q2 .* Fpu;
86 Pv = p1 .* Ftv + p2 .* Fpv;
87 Qv = p2 .* Ftv + q2 .* Fpv;
a4fee9cBatch independent transforms through one Legendre dispatchDan Fortunato 88 [PAu, QAu, PAv, QAv] = analys(Pu, Qu, Pv, Qv);
89 Pcu = PAu .* filt;
90 Qcu = QAu .* filt;
91 Pcv = PAv .* filt;
92 Qcv = QAv .* filt;
93 scu = dthetac(Pcu) + dphic(Qcu);
591a4f5Reduce the Laplace-Beltrami matvec to 6 transforms per species per iterationDan Fortunato 94 scv = dthetac(Pcv) + dphic(Qcv);
a4fee9cBatch independent transforms through one Legendre dispatchDan Fortunato 95 [Lu, Lv] = synth(scu, scv);
96 lapu = r .* Lu;
97 lapv = r .* Lv;
98 [LAu, LAv] = analys(lapu, lapv);
e4d6a3bPrecondition with the operator's symbol; project the correction onto the bandDan Fortunato 99 dLu = (LAu + lamJ .* Un) .* filt;
100 dLv = (LAv + lamJ .* Vn) .* filt;
8144287WIP: added code for Laplace-Beltrami operator evaluation on smooth genus-0 surfaceOwen Melia 101
e4d6a3bPrecondition with the operator's symbol; project the correction onto the bandDan Fortunato 102 Un = (Bu + (dt * D1) * dLu) ./ (1 + (dt * D1) * lamJ);
103 Vn = (Bv + (dt * D2) * dLv) ./ (1 + (dt * D2) * lamJ);
105end