62e6cc4Simplify UI text and built-in .m script commentsJeremy Magland 1% Brusselator reaction-diffusion on a closed surface.
3% du/dt = D1*lap_g(u) + A - (B+1)*u + u^2*v
4% dv/dt = D2*lap_g(v) + B*u - u^2*v
5%
a4fee9cBatch independent transforms through one Legendre dispatchDan Fortunato 6% Same scheme as models/schnakenberg.m, including the grouped transforms:
3d078cfSplit the flux-form divergence against the round sphereDan Fortunato 7% [a, b] = synth(x, y) runs the group as batched Legendre dispatches, and the
8% sphere-split flux divergence that keeps r ~ 1/sin^2(theta) off the round
9% sphere's share of the operator.
0ae15cfSeed runs from smooth random fields, and add the blob geometryDan Fortunato 11% Seeded from a smooth random field -- see models/schnakenberg.m.
12function [U, V, u, v] = init(lam3, gx, gy, gz, A, B)
13 f = randnfun3(lam3, gx, gy, gz);
14 [U, V] = analys(A + 0.01*f, (B / A) * ones(numel(f), 1));
a4fee9cBatch independent transforms through one Legendre dispatchDan Fortunato 15 [u, v] = synth(U, V);
3d078cfSplit the flux-form divergence against the round sphereDan Fortunato 18function [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 19 [u, v] = synth(U, V);
a4fee9cBatch independent transforms through one Legendre dispatchDan Fortunato 22 ru = A - (B + 1) * u + uuv;
23 rv = B * u - uuv;
24 [Ru, Rv] = analys(ru, rv);
25 Bu = U + dt * Ru;
26 Bv = V + dt * Rv;
e4d6a3bPrecondition with the operator's symbol; project the correction onto the bandDan Fortunato 28 % Mean-J preconditioning -- see models/schnakenberg.m.
29 lamJ = lam ./ jhat;
30 Un = Bu ./ (1 + (dt * D1) * lamJ);
31 Vn = Bv ./ (1 + (dt * D2) * lamJ);
33 for k = 1:niter
591a4f5Reduce the Laplace-Beltrami matvec to 6 transforms per species per iterationDan Fortunato 34 % dlap = lap_g - lap_s, evaluated at the current iterate in flux form
35 % (see models/schnakenberg.m, docs/richardson-iteration.md and
a4fee9cBatch independent transforms through one Legendre dispatchDan Fortunato 36 % docs/reduced-transforms.md for the derivation and the ordering).
8144287WIP: added code for Laplace-Beltrami operator evaluation on smooth genus-0 surfaceOwen Melia 37 Fu = Un .* filt;
39 vtu = dthetac(Fu);
40 vpu = dphic(Fu);
41 vtv = dthetac(Fv);
42 vpv = dphic(Fv);
3d078cfSplit the flux-form divergence against the round sphereDan Fortunato 43 [Ftu, Fpu, Ftv, Fpv, Su, Sv] = synth(vtu, vpu, vtv, vpv, lam .* Fu, lam .* Fv);
44 Pu = dp1 .* Ftu + p2 .* Fpu;
45 Qu = p2 .* Ftu + dq2 .* Fpu;
46 Pv = dp1 .* Ftv + p2 .* Fpv;
47 Qv = p2 .* Ftv + dq2 .* Fpv;
a4fee9cBatch independent transforms through one Legendre dispatchDan Fortunato 49 Pcu = PAu .* filt;
50 Pcv = PAv .* filt;
52 scv = dthetac(Pcv);
a4fee9cBatch independent transforms through one Legendre dispatchDan Fortunato 53 [Lu, Lv] = synth(scu, scv);
55 dQv = dphig(Qv);
3d078cfSplit the flux-form divergence against the round sphereDan Fortunato 56 lapu = r .* (Lu + dQu) - jinv .* Su;
57 lapv = r .* (Lv + dQv) - jinv .* Sv;
a4fee9cBatch independent transforms through one Legendre dispatchDan Fortunato 58 [LAu, LAv] = analys(lapu, lapv);
e4d6a3bPrecondition with the operator's symbol; project the correction onto the bandDan Fortunato 59 dLu = (LAu + lamJ .* Un) .* filt;
60 dLv = (LAv + lamJ .* Vn) .* filt;
8144287WIP: added code for Laplace-Beltrami operator evaluation on smooth genus-0 surfaceOwen Melia 61
e4d6a3bPrecondition with the operator's symbol; project the correction onto the bandDan Fortunato 62 Un = (Bu + (dt * D1) * dLu) ./ (1 + (dt * D1) * lamJ);
63 Vn = (Bv + (dt * D2) * dLv) ./ (1 + (dt * D2) * lamJ);
65end