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.
0af3386Reaction-diffusion on spherical-harmonic surfacesJeremy Magland 11
12function [U, V, u, v] = init(noise, a, b)
13 us = a + b;
14 vs = b / (us * us);
15 U = analys(us + noise);
16 V = analys(vs * ones(numel(noise), 1));
17 u = synth(U);
18 v = synth(V);
19end
8144287WIP: added code for Laplace-Beltrami operator evaluation on smooth genus-0 surfaceOwen Melia 21function [Un, Vn, u, v] = step(U, V, lam, filt, gx, gy, gz, Vtx, Vty, Vtz, Vpx, Vpy, Vpz, a, b, D1, D2, dt, niter)
23 v = synth(V);
24 uuv = u .* u .* v;
62e6cc4Simplify UI text and built-in .m script commentsJeremy Magland 26 % Right-hand side of the implicit solve (I - dt*D*lap_g) Unew = B.
0af3386Reaction-diffusion on spherical-harmonic surfacesJeremy Magland 27 Bu = U + dt * analys(a - u + uuv);
28 Bv = V + dt * analys(b - uuv);
62e6cc4Simplify UI text and built-in .m script commentsJeremy Magland 30 % Round-sphere solve, then iterate the geometric correction.
0af3386Reaction-diffusion on spherical-harmonic surfacesJeremy Magland 31 Un = Bu ./ (1 + (dt * D1) * lam);
32 Vn = Bv ./ (1 + (dt * D2) * lam);
34 for k = 1:niter
8144287WIP: added code for Laplace-Beltrami operator evaluation on smooth genus-0 surfaceOwen Melia 35 % dlap = lap_g - lap_s, evaluated at the current iterate (Algorithm 3 of
36 % evolving_surface/notes/algos.tex): surface gradient of the field,
37 % contracted through the inverse metric quantities Vt*/Vp*; each
38 % Cartesian component re-analysed and differentiated again; recombined
39 % into the surface divergence. lam.*Un adds back -lap_s(Un), since lam
40 % holds +l(l+1). filt zeroes the top two degrees, where the theta/phi
41 % derivative recurrences cannot exactly represent a derivative. See
42 % docs/richardson-iteration.md.
43 Fu = Un .* filt;
44 Ftu = dtheta(Fu);
45 Fpu = dphi(Fu);
46 dux = Ftu .* Vtx + Fpu .* Vpx;
47 duy = Ftu .* Vty + Fpu .* Vpy;
48 duz = Ftu .* Vtz + Fpu .* Vpz;
49 cux = analys(dux) .* filt;
50 cuy = analys(duy) .* filt;
51 cuz = analys(duz) .* filt;
52 Ftcux = dtheta(cux);
53 Fpcux = dphi(cux);
54 Ftcuy = dtheta(cuy);
55 Fpcuy = dphi(cuy);
56 Ftcuz = dtheta(cuz);
57 Fpcuz = dphi(cuz);
58 lapu = Ftcux .* Vtx + Fpcux .* Vpx;
59 lapu = lapu + Ftcuy .* Vty;
60 lapu = lapu + Fpcuy .* Vpy;
61 lapu = lapu + Ftcuz .* Vtz;
62 lapu = lapu + Fpcuz .* Vpz;
63 dLu = analys(lapu) + lam .* Un;
65 Fv = Vn .* filt;
66 Ftv = dtheta(Fv);
67 Fpv = dphi(Fv);
68 dvx = Ftv .* Vtx + Fpv .* Vpx;
69 dvy = Ftv .* Vty + Fpv .* Vpy;
70 dvz = Ftv .* Vtz + Fpv .* Vpz;
71 cvx = analys(dvx) .* filt;
72 cvy = analys(dvy) .* filt;
73 cvz = analys(dvz) .* filt;
74 Ftcvx = dtheta(cvx);
75 Fpcvx = dphi(cvx);
76 Ftcvy = dtheta(cvy);
77 Fpcvy = dphi(cvy);
78 Ftcvz = dtheta(cvz);
79 Fpcvz = dphi(cvz);
80 lapv = Ftcvx .* Vtx + Fpcvx .* Vpx;
81 lapv = lapv + Ftcvy .* Vty;
82 lapv = lapv + Fpcvy .* Vpy;
83 lapv = lapv + Ftcvz .* Vtz;
84 lapv = lapv + Fpcvz .* Vpz;
85 dLv = analys(lapv) + lam .* Vn;
0af3386Reaction-diffusion on spherical-harmonic surfacesJeremy Magland 87 Un = (Bu + (dt * D1) * dLu) ./ (1 + (dt * D1) * lam);
88 Vn = (Bv + (dt * D2) * dLv) ./ (1 + (dt * D2) * lam);
89 end
90end