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 step forms the
7% right-hand side B of the linear system (I - dt*D*lap_g) Unew = B, and hands
8% the solve to solve(...) — the solver the app's selector picks (richardson,
9% bicgstab or gmres, see solvers/), every one applying the operator's
10% geometric part through lib/dlap.m. A model may also name a solver directly
11% in these two call lines. Grid fields are npts x 1; spectral fields are
12% real 2 x nlm. See docs/richardson-iteration.md.
14function [U, V, u, v] = init(noise, a, b)
15 us = a + b;
16 vs = b / (us * us);
17 U = analys(us + noise);
18 V = analys(vs * ones(numel(noise), 1));
19 u = synth(U);
20 v = synth(V);
21end
23function [Un, Vn, u, v] = step(U, V, lam, filt, wlm, gx, gy, gz, Vtx, Vty, Vtz, Vpx, Vpy, Vpz, a, b, D1, D2, dt, nlm, niter)
24 u = synth(U);
25 v = synth(V);
26 uuv = u .* u .* v;
28 % Right-hand side of the implicit solve (I - dt*D*lap_g) Unew = B.
29 Bu = U + dt * analys(a - u + uuv);
30 Bv = V + dt * analys(b - uuv);
32 % The species diffuse independently, so each gets its own solve.
33 Un = solve(Bu, dt * D1, lam, filt, wlm, Vtx, Vty, Vtz, Vpx, Vpy, Vpz, nlm, niter);
34 Vn = solve(Bv, dt * D2, lam, filt, wlm, Vtx, Vty, Vtz, Vpx, Vpy, Vpz, nlm, niter);
35end