1% The geometric part of the surface Laplace-Beltrami operator:
2%
3% dlap = lap_g - lap_s
4%
5% applied to a spectral field F, where lap_g is the surface's operator and
6% lap_s the round sphere's (diagonal, eigenvalues -l(l+1)). This is the piece
7% of the implicit solve (I - dt*D*lap_g) X = B that a solver re-evaluates
8% each iteration — the round-sphere part it inverts exactly. The operator is
9% linear in F; the surface enters only through the inverse metric quantities
10% Vt*/Vp* (src/geom/metric.ts), so swapping the geometry changes no code.
11%
12% Algorithm 3 of evolving_surface/notes/algos.tex: surface gradient of the
13% field (theta/phi derivatives contracted through Vt*/Vp*); each Cartesian
14% component re-analysed and differentiated again; recombined into the surface
15% divergence. lam .* F adds back -lap_s(F), since lam holds +l(l+1). filt
16% zeroes the top two degrees, where the theta/phi derivative recurrences
17% cannot exactly represent a derivative. See docs/richardson-iteration.md.
18%
19% Spectral fields are real 2 x nlm; the intermediate d*/L fields live on the
20% npts x 1 grid.
22function dL = dlap(F, filt, Vtx, Vty, Vtz, Vpx, Vpy, Vpz, lam)
23 G = F .* filt;
24 Ft = dtheta(G);
25 Fp = dphi(G);
26 dx = Ft .* Vtx + Fp .* Vpx;
27 dy = Ft .* Vty + Fp .* Vpy;
28 dz = Ft .* Vtz + Fp .* Vpz;
29 cx = analys(dx) .* filt;
30 cy = analys(dy) .* filt;
31 cz = analys(dz) .* filt;
32 Ftcx = dtheta(cx);
33 Fpcx = dphi(cx);
34 Ftcy = dtheta(cy);
35 Fpcy = dphi(cy);
36 Ftcz = dtheta(cz);
37 Fpcz = dphi(cz);
38 L = Ftcx .* Vtx + Fpcx .* Vpx;
39 L = L + Ftcy .* Vty;
40 L = L + Fpcy .* Vpy;
41 L = L + Ftcz .* Vtz;
42 L = L + Fpcz .* Vpz;
43 dL = analys(L) + lam .* F;
44end