591a4f5Reduce the Laplace-Beltrami matvec to 6 transforms per species per iterationDan Fortunato 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%
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;
10% spectral fields are real 2 x nlm. See docs/richardson-iteration.md.
11%
12% This is the 12-transform reference: the correction evaluates lap_g in
13% Cartesian-gradient form (Algorithm 4 of evolving_surface/notes/algos.tex),
14% carrying grad_g(u) as three ambient components through the inverse metric
15% quantities Vt*/Vp*. models/schnakenberg.m computes the same operator in
16% flux form with 6 transforms per species per iteration
17% (docs/reduced-transforms.md); this variant is kept live
18% for A/B comparison, in the app and in the tests.
0ae15cfSeed runs from smooth random fields, and add the blob geometryDan Fortunato 20% Seeded from a smooth random field -- see models/schnakenberg.m.
21function [U, V, u, v] = init(lam3, gx, gy, gz, a, b)
22 f = randnfun3(lam3, gx, gy, gz);
591a4f5Reduce the Laplace-Beltrami matvec to 6 transforms per species per iterationDan Fortunato 23 us = a + b;
24 vs = b / (us * us);
0ae15cfSeed runs from smooth random fields, and add the blob geometryDan Fortunato 25 U = analys(us + 0.01*f);
26 V = analys(vs * ones(numel(f), 1));
591a4f5Reduce the Laplace-Beltrami matvec to 6 transforms per species per iterationDan Fortunato 27 u = synth(U);
28 v = synth(V);
29end
e4d6a3bPrecondition with the operator's symbol; project the correction onto the bandDan Fortunato 31function [Un, Vn, u, v] = step(U, V, lam, filt, gx, gy, gz, Vtx, Vty, Vtz, Vpx, Vpy, Vpz, jhat, a, b, D1, D2, dt, niter)
591a4f5Reduce the Laplace-Beltrami matvec to 6 transforms per species per iterationDan Fortunato 32 u = synth(U);
33 v = synth(V);
34 uuv = u .* u .* v;
36 % Right-hand side of the implicit solve (I - dt*D*lap_g) Unew = B.
37 Bu = U + dt * analys(a - u + uuv);
38 Bv = V + dt * analys(b - uuv);
e4d6a3bPrecondition with the operator's symbol; project the correction onto the bandDan Fortunato 40 % Mean-J preconditioned solve (see models/schnakenberg.m), then iterate
41 % the geometric correction.
42 lamJ = lam ./ jhat;
43 Un = Bu ./ (1 + (dt * D1) * lamJ);
44 Vn = Bv ./ (1 + (dt * D2) * lamJ);
46 for k = 1:niter
47 % dlap = lap_g - lap_s, evaluated at the current iterate (Algorithm 3 of
48 % evolving_surface/notes/algos.tex): surface gradient of the field,
49 % contracted through the inverse metric quantities Vt*/Vp*; each
50 % Cartesian component re-analysed and differentiated again; recombined
51 % into the surface divergence. lam.*Un adds back -lap_s(Un), since lam
52 % holds +l(l+1). filt zeroes the top two degrees, where the theta/phi
53 % derivative recurrences cannot exactly represent a derivative. See
54 % docs/richardson-iteration.md.
55 Fu = Un .* filt;
56 Ftu = dtheta(Fu);
57 Fpu = dphi(Fu);
58 dux = Ftu .* Vtx + Fpu .* Vpx;
59 duy = Ftu .* Vty + Fpu .* Vpy;
60 duz = Ftu .* Vtz + Fpu .* Vpz;
61 cux = analys(dux) .* filt;
62 cuy = analys(duy) .* filt;
63 cuz = analys(duz) .* filt;
64 Ftcux = dtheta(cux);
65 Fpcux = dphi(cux);
66 Ftcuy = dtheta(cuy);
67 Fpcuy = dphi(cuy);
68 Ftcuz = dtheta(cuz);
69 Fpcuz = dphi(cuz);
70 lapu = Ftcux .* Vtx + Fpcux .* Vpx;
71 lapu = lapu + Ftcuy .* Vty;
72 lapu = lapu + Fpcuy .* Vpy;
73 lapu = lapu + Ftcuz .* Vtz;
74 lapu = lapu + Fpcuz .* Vpz;
e4d6a3bPrecondition with the operator's symbol; project the correction onto the bandDan Fortunato 75 dLu = (analys(lapu) + lamJ .* Un) .* filt;
77 Fv = Vn .* filt;
78 Ftv = dtheta(Fv);
79 Fpv = dphi(Fv);
80 dvx = Ftv .* Vtx + Fpv .* Vpx;
81 dvy = Ftv .* Vty + Fpv .* Vpy;
82 dvz = Ftv .* Vtz + Fpv .* Vpz;
83 cvx = analys(dvx) .* filt;
84 cvy = analys(dvy) .* filt;
85 cvz = analys(dvz) .* filt;
86 Ftcvx = dtheta(cvx);
87 Fpcvx = dphi(cvx);
88 Ftcvy = dtheta(cvy);
89 Fpcvy = dphi(cvy);
90 Ftcvz = dtheta(cvz);
91 Fpcvz = dphi(cvz);
92 lapv = Ftcvx .* Vtx + Fpcvx .* Vpx;
93 lapv = lapv + Ftcvy .* Vty;
94 lapv = lapv + Fpcvy .* Vpy;
95 lapv = lapv + Ftcvz .* Vtz;
96 lapv = lapv + Fpcvz .* Vpz;
e4d6a3bPrecondition with the operator's symbol; project the correction onto the bandDan Fortunato 97 dLv = (analys(lapv) + lamJ .* Vn) .* filt;
e4d6a3bPrecondition with the operator's symbol; project the correction onto the bandDan Fortunato 99 Un = (Bu + (dt * D1) * dLu) ./ (1 + (dt * D1) * lamJ);
100 Vn = (Bv + (dt * D2) * dLv) ./ (1 + (dt * D2) * lamJ);
591a4f5Reduce the Laplace-Beltrami matvec to 6 transforms per species per iterationDan Fortunato 101 end
102end