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%
3d078cfSplit the flux-form divergence against the round sphereDan Fortunato 12% The correction evaluates lap_g in flux form -- 7 transforms per species
591a4f5Reduce the Laplace-Beltrami matvec to 6 transforms per species per iterationDan Fortunato 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.
18% The flux divergence is split against the round sphere: the sphere's share
19% of it is -jinv*lap_s(u), exact in spectral space, and only the geometry
20% *deviation* meets r ~ 1/sin^2(theta). Without that split the concentrated
21% division amplifies the polar roundoff of the whole flux, and since a
22% Turing pattern is seeded by whatever is largest in its unstable band, the
23% amplified polar noise -- static, and re-injected every step -- picks the
24% nucleation site and grows a spot at the pole. See docs/reduced-transforms.md
25% Sec 5.
0ae15cfSeed runs from smooth random fields, and add the blob geometryDan Fortunato 27% The uniform steady state, perturbed by a smooth random field: chebfun's
28% randnfun3 on the surface's bounding box, restricted to the surface by
29% evaluating it at the grid points -- the way surfacefun seeds a run. lam3
30% is its wavelength; the draw is seeded on the host, the sum over its
31% Fourier modes runs on the GPU (src/mgpu/randnfun3.ts).
32function [U, V, u, v] = init(lam3, gx, gy, gz, a, b)
33 f = randnfun3(lam3, gx, gy, gz);
35 vs = b / (us * us);
0ae15cfSeed runs from smooth random fields, and add the blob geometryDan Fortunato 36 [U, V] = analys(us + 0.01*f, vs * ones(numel(f), 1));
a4fee9cBatch independent transforms through one Legendre dispatchDan Fortunato 37 [u, v] = synth(U, V);
3d078cfSplit the flux-form divergence against the round sphereDan Fortunato 40function [Un, Vn, u, v] = step(U, V, lam, filt, gx, gy, gz, p2, r, dp1, dq2, jinv, jhat, a, b, D1, D2, dt, niter)
a4fee9cBatch independent transforms through one Legendre dispatchDan Fortunato 41 % Grouped transforms -- [a, b] = synth(x, y) -- are explicit batching:
42 % output k is the transform of input k, and the whole group runs as one
43 % batched Legendre dispatch, or as many as the device's lane width allows
44 % (src/mgpu/plan.ts, materializeTransforms). The grouping is a promise of
45 % independence, never of a lane width, so the same source runs anywhere.
46 [u, v] = synth(U, V);
62e6cc4Simplify UI text and built-in .m script commentsJeremy Magland 49 % Right-hand side of the implicit solve (I - dt*D*lap_g) Unew = B.
51 rv = b - uuv;
52 [Ru, Rv] = analys(ru, rv);
53 Bu = U + dt * Ru;
54 Bv = V + dt * Rv;
0af3386Reaction-diffusion on spherical-harmonic surfacesJeremy Magland 55
e4d6a3bPrecondition with the operator's symbol; project the correction onto the bandDan Fortunato 56 % Preconditioned solve, then iterate the geometric correction. jhat is
57 % the host's minimax scale over the operator's symbol eigenvalues mu(x)
58 % -- the inverse squared principal stretches of the embedding, direction
59 % included (src/geom/geometry.ts, Jhat): preconditioning with lam/jhat
60 % contracts every mode and direction at rate
61 % (muMax - muMin)/(muMax + muMin) < 1 on any surface, where the plain
62 % lam (jhat = 1) diverges wherever mu > 2 -- docs/reduced-transforms.md
63 % Sec 10. The answer never depends on jhat (the lamJ term added inside
64 % dLu is the term divided back out); only the convergence rate does. On
65 % the sphere mu = 1 and lamJ = lam.
66 lamJ = lam ./ jhat;
67 Un = Bu ./ (1 + (dt * D1) * lamJ);
68 Vn = Bv ./ (1 + (dt * D2) * lamJ);
0af3386Reaction-diffusion on spherical-harmonic surfacesJeremy Magland 69
70 for k = 1:niter
591a4f5Reduce the Laplace-Beltrami matvec to 6 transforms per species per iterationDan Fortunato 71 % dlap = lap_g - lap_s at the current iterate, in flux form
a4fee9cBatch independent transforms through one Legendre dispatchDan Fortunato 72 % (docs/reduced-transforms.md Sec 4). The sin-weighted derivatives
3d078cfSplit the flux-form divergence against the round sphereDan Fortunato 73 % A = sin(theta)*dtheta(u) and B = dphi(u) -- both smooth on the sphere,
a4fee9cBatch independent transforms through one Legendre dispatchDan Fortunato 74 % synthesized straight from the dthetac/dphic coefficient shuffles --
3d078cfSplit the flux-form divergence against the round sphereDan Fortunato 75 % are combined pointwise through the precomputed weights into two
76 % fluxes P,Q, also smooth. The theta flux P goes back to
0d99c91Differentiate the phi flux in grid spaceDan Fortunato 77 % coefficients, through the same shuffle again, and is synthesized as
78 % sin(theta)*dtheta(P); the phi flux Q never leaves the grid -- d/dphi
79 % is diagonal in the Fourier index, so dphig differentiates it with two
80 % FFT stages and no Legendre work (masking m past filt's reach). Their
81 % sum, scaled by r, is lap_g(u). The only division by sin(theta)
3d078cfSplit the flux-form divergence against the round sphereDan Fortunato 82 % anywhere is folded into the weights at precompute time.
83 %
84 % The weights here are the *sphere-subtracted* ones: p1 = 1 + dp1 and
85 % q2 = 1 + dq2 (p2 is zero on the sphere already), so P,Q below are the
86 % deviation fluxes P' = P - A, Q' = Q - B. What that leaves out is the
87 % round sphere's own divergence, sin(theta)*dtheta(A) + dphi(B) =
88 % -sin^2(theta)*lap_s(u), which needs no flux machinery at all: lap_s is
89 % diagonal, so it is -lam.*Fu synthesized once (S below, riding along in
90 % the gradient's batched synthesis) and scaled by the bounded
91 % jinv = 1/J = r*sin^2(theta). r therefore multiplies only the deviation
92 % -- the difference between this and multiplying the whole flux is two
93 % orders of magnitude of polar roundoff, and it is what keeps a pattern
94 % from nucleating at the pole (src/geom/geometry.ts, dp1/dq2/jinv).
e4d6a3bPrecondition with the operator's symbol; project the correction onto the bandDan Fortunato 95 % lamJ.*Un adds back the preconditioner's -lap_s(Un)/jhat, since lam
96 % holds +l(l+1). filt zeroes the top two degrees, where the derivative
97 % recurrences cannot exactly represent a derivative -- and the correction
98 % itself is projected onto the same band (algos.tex Algorithm 5 zeroes
99 % the same coefficients): without that, each iteration replaces a bit
100 % more of the top degrees' implicit diffusion with nothing (their fixed
101 % point is the undiffused Bu), and the two species un-diffuse at
102 % different rates -- a spurious Turing band at the band edge.
3d078cfSplit the flux-form divergence against the round sphereDan Fortunato 104 % The two species share each grouped call: the six gradient-and-sphere
0d99c91Differentiate the phi flux in grid spaceDan Fortunato 105 % syntheses, the two theta-flux analyses, the two divergence syntheses
106 % and the two final analyses each run as one batched dispatch.
8144287WIP: added code for Laplace-Beltrami operator evaluation on smooth genus-0 surfaceOwen Melia 107 Fu = Un .* filt;
109 vtu = dthetac(Fu);
110 vpu = dphic(Fu);
111 vtv = dthetac(Fv);
112 vpv = dphic(Fv);
3d078cfSplit the flux-form divergence against the round sphereDan Fortunato 113 [Ftu, Fpu, Ftv, Fpv, Su, Sv] = synth(vtu, vpu, vtv, vpv, lam .* Fu, lam .* Fv);
114 Pu = dp1 .* Ftu + p2 .* Fpu;
115 Qu = p2 .* Ftu + dq2 .* Fpu;
116 Pv = dp1 .* Ftv + p2 .* Fpv;
117 Qv = p2 .* Ftv + dq2 .* Fpv;
a4fee9cBatch independent transforms through one Legendre dispatchDan Fortunato 119 Pcu = PAu .* filt;
120 Pcv = PAv .* filt;
122 scv = dthetac(Pcv);
a4fee9cBatch independent transforms through one Legendre dispatchDan Fortunato 123 [Lu, Lv] = synth(scu, scv);
125 dQv = dphig(Qv);
3d078cfSplit the flux-form divergence against the round sphereDan Fortunato 126 lapu = r .* (Lu + dQu) - jinv .* Su;
127 lapv = r .* (Lv + dQv) - jinv .* Sv;
a4fee9cBatch independent transforms through one Legendre dispatchDan Fortunato 128 [LAu, LAv] = analys(lapu, lapv);
e4d6a3bPrecondition with the operator's symbol; project the correction onto the bandDan Fortunato 129 dLu = (LAu + lamJ .* Un) .* filt;
130 dLv = (LAv + lamJ .* Vn) .* filt;
8144287WIP: added code for Laplace-Beltrami operator evaluation on smooth genus-0 surfaceOwen Melia 131
e4d6a3bPrecondition with the operator's symbol; project the correction onto the bandDan Fortunato 132 Un = (Bu + (dt * D1) * dLu) ./ (1 + (dt * D1) * lamJ);
133 Vn = (Bv + (dt * D2) * dLv) ./ (1 + (dt * D2) * lamJ);
135end