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);
21 U = analys(us + noise);
22 V = analys(vs * ones(numel(noise), 1));
23 u = synth(U);
24 v = synth(V);
25end
591a4f5Reduce the Laplace-Beltrami matvec to 6 transforms per species per iterationDan Fortunato 27function [Un, Vn, u, v] = step(U, V, lam, filt, gx, gy, gz, p1, p2, q2, r, a, b, D1, D2, dt, niter)
29 v = synth(V);
30 uuv = u .* u .* v;
62e6cc4Simplify UI text and built-in .m script commentsJeremy Magland 32 % Right-hand side of the implicit solve (I - dt*D*lap_g) Unew = B.
0af3386Reaction-diffusion on spherical-harmonic surfacesJeremy Magland 33 Bu = U + dt * analys(a - u + uuv);
34 Bv = V + dt * analys(b - uuv);
62e6cc4Simplify UI text and built-in .m script commentsJeremy Magland 36 % Round-sphere solve, then iterate the geometric correction.
0af3386Reaction-diffusion on spherical-harmonic surfacesJeremy Magland 37 Un = Bu ./ (1 + (dt * D1) * lam);
38 Vn = Bv ./ (1 + (dt * D2) * lam);
40 for k = 1:niter
591a4f5Reduce the Laplace-Beltrami matvec to 6 transforms per species per iterationDan Fortunato 41 % dlap = lap_g - lap_s at the current iterate, in flux form
42 % (docs/reduced-transforms.md Sec 4). The sin-weighted
43 % derivatives sin(theta)*dtheta(u) and dphi(u) -- both smooth on the
44 % sphere, synthesized straight from the dthetac/dphic coefficient
45 % shuffles -- are combined pointwise through the precomputed weights
46 % p1,p2,q2 into two fluxes P,Q, also smooth. Their coefficients are then
47 % pushed through the *same* shuffles again and summed before the one
48 % synthesis of the divergence, which r scales into lap_g(u). The only
49 % division by sin(theta) anywhere is folded into p1,p2,q2,r at precompute
50 % time. lam.*Un adds back -lap_s(Un), since lam holds +l(l+1). filt
51 % zeroes the top two degrees, where the derivative recurrences cannot
52 % exactly represent a derivative.
8144287WIP: added code for Laplace-Beltrami operator evaluation on smooth genus-0 surfaceOwen Melia 53 Fu = Un .* filt;
591a4f5Reduce the Laplace-Beltrami matvec to 6 transforms per species per iterationDan Fortunato 54 Ftu = synth(dthetac(Fu));
55 Fpu = synth(dphic(Fu));
56 Pu = p1 .* Ftu + p2 .* Fpu;
57 Qu = p2 .* Ftu + q2 .* Fpu;
58 Pcu = analys(Pu) .* filt;
59 Qcu = analys(Qu) .* filt;
60 scu = dthetac(Pcu) + dphic(Qcu);
61 lapu = r .* synth(scu);
8144287WIP: added code for Laplace-Beltrami operator evaluation on smooth genus-0 surfaceOwen Melia 62 dLu = analys(lapu) + lam .* Un;
64 Fv = Vn .* filt;
591a4f5Reduce the Laplace-Beltrami matvec to 6 transforms per species per iterationDan Fortunato 65 Ftv = synth(dthetac(Fv));
66 Fpv = synth(dphic(Fv));
67 Pv = p1 .* Ftv + p2 .* Fpv;
68 Qv = p2 .* Ftv + q2 .* Fpv;
69 Pcv = analys(Pv) .* filt;
70 Qcv = analys(Qv) .* filt;
71 scv = dthetac(Pcv) + dphic(Qcv);
72 lapv = r .* synth(scv);
8144287WIP: added code for Laplace-Beltrami operator evaluation on smooth genus-0 surfaceOwen Melia 73 dLv = analys(lapv) + lam .* Vn;
0af3386Reaction-diffusion on spherical-harmonic surfacesJeremy Magland 75 Un = (Bu + (dt * D1) * dLu) ./ (1 + (dt * D1) * lam);
76 Vn = (Bv + (dt * D2) * dLv) ./ (1 + (dt * D2) * lam);
77 end
78end