2dedc35turing-sphere: reaction-diffusion on the sphere, spectral solver on WebGPUJeremy Magland 1/**
2 * Double-precision reference implementation of the scalar spherical
3 * harmonic transform, by direct summation. Slow (O(nlat*nlm) Legendre +
4 * O(nlat*nphi*mmax) Fourier) but simple, and serves as ground truth for
5 * validating the fp32 WebGPU implementation.
6 *
7 * Conventions are identical to the GPU path (see layout.ts).
8 */
9import { gaussNodesWeights } from './gauss.ts';
10import { legendreCoeffs, legendreRow, type LegendreCoeffs } from './coeffs.ts';
11import { lmIndex, nlmCalc, validateConfig, type ShtConfig } from './layout.ts';
13export class ShtReference {
14 readonly cfg: ShtConfig;
15 readonly nlm: number;
16 readonly ct: Float64Array;
17 readonly st: Float64Array;
18 readonly wg: Float64Array; // Gauss weights (for integral over cos(theta))
19 readonly coeffs: LegendreCoeffs;
21 constructor(cfg: ShtConfig) {
22 validateConfig(cfg);
23 this.cfg = cfg;
24 this.nlm = nlmCalc(cfg.lmax, cfg.mmax);
25 const { x, w } = gaussNodesWeights(cfg.nlat);
26 this.ct = x;
27 this.wg = w;
28 this.st = new Float64Array(cfg.nlat);
29 for (let i = 0; i < cfg.nlat; i++) this.st[i] = Math.sqrt(1 - x[i] * x[i]);
30 this.coeffs = legendreCoeffs(cfg.lmax, cfg.mmax);
31 }
33 /**
34 * Legendre stage of the synthesis: F_m(theta_i) = sum_l Q_lm ytilde_l^m(theta_i).
35 * Returns complex array indexed [m * nlat + ilat], interleaved re/im.
36 */
37 legendreSynth(qlm: ArrayLike<number>): Float64Array {
38 const { lmax, mmax, nlat } = this.cfg;
39 const fm = new Float64Array(2 * (mmax + 1) * nlat);
40 const row = new Float64Array(lmax + 1);
41 for (let i = 0; i < nlat; i++) {
42 for (let m = 0; m <= mmax; m++) {
43 legendreRow(this.coeffs, lmax, m, this.ct[i], this.st[i], row);
44 let re = 0, im = 0;
45 const base = lmIndex(lmax, m, m);
46 for (let l = m; l <= lmax; l++) {
47 const y = row[l - m];
48 re += y * qlm[2 * (base + l - m)];
49 im += y * qlm[2 * (base + l - m) + 1];
50 }
51 const o = 2 * (m * nlat + i);
52 fm[o] = re;
53 fm[o + 1] = im;
54 }
55 }
56 return fm;
57 }
59 /** Full synthesis: spectral -> spatial grid [ilat * nphi + iphi]. */
60 synth(qlm: ArrayLike<number>): Float64Array {
61 const { mmax, nlat, nphi } = this.cfg;
62 const fm = this.legendreSynth(qlm);
63 const spat = new Float64Array(nlat * nphi);
64 for (let i = 0; i < nlat; i++) {
65 for (let j = 0; j < nphi; j++) {
66 const phi = (2 * Math.PI * j) / nphi;
67 let v = fm[2 * (0 * nlat + i)]; // m=0: real part (imag must be 0)
68 for (let m = 1; m <= mmax; m++) {
69 const o = 2 * (m * nlat + i);
70 const c = Math.cos(m * phi);
71 const s = Math.sin(m * phi);
72 v += 2 * (fm[o] * c - fm[o + 1] * s);
73 }
74 spat[i * nphi + j] = v;
75 }
76 }
77 return spat;
78 }
80 /** Full analysis: spatial grid -> spectral coefficients (interleaved re/im). */
81 analys(spat: ArrayLike<number>): Float64Array {
82 const { lmax, mmax, nlat, nphi } = this.cfg;
83 const qlm = new Float64Array(2 * this.nlm);
84 const row = new Float64Array(lmax + 1);
85 // forward Fourier: G_m(theta_i) = (2*pi/nphi) * sum_j f_ij e^{-i m phi_j}
86 const gm = new Float64Array(2 * (mmax + 1) * nlat);
87 for (let i = 0; i < nlat; i++) {
88 for (let m = 0; m <= mmax; m++) {
89 let re = 0, im = 0;
90 for (let j = 0; j < nphi; j++) {
91 const phi = (2 * Math.PI * j) / nphi;
92 const f = spat[i * nphi + j];
93 re += f * Math.cos(m * phi);
94 im -= f * Math.sin(m * phi);
95 }
96 const o = 2 * (m * nlat + i);
97 const norm = (2 * Math.PI) / nphi;
98 gm[o] = re * norm;
99 gm[o + 1] = im * norm;
100 }
101 }
102 // Legendre stage with Gauss quadrature: Q_lm = sum_i w_i ytilde_l^m(theta_i) G_m(theta_i)
103 for (let m = 0; m <= mmax; m++) {
104 const base = lmIndex(lmax, m, m);
105 for (let i = 0; i < nlat; i++) {
106 legendreRow(this.coeffs, lmax, m, this.ct[i], this.st[i], row);
107 const o = 2 * (m * nlat + i);
108 const wr = this.wg[i] * gm[o];
109 const wi = this.wg[i] * gm[o + 1];
110 for (let l = m; l <= lmax; l++) {
111 const y = row[l - m];
112 qlm[2 * (base + l - m)] += y * wr;
113 qlm[2 * (base + l - m) + 1] += y * wi;
114 }
115 }
116 }
117 return qlm;
118 }
119}
121/** Random band-limited spectrum for testing (m=0 imaginary parts zeroed). */
122export function randomSpectrum(cfg: ShtConfig, seed = 12345): Float32Array {
123 const nlm = nlmCalc(cfg.lmax, cfg.mmax);
124 const q = new Float32Array(2 * nlm);
125 let s = seed >>> 0;
126 const rnd = () => {
127 // xorshift32
128 s ^= s << 13; s >>>= 0;
129 s ^= s >> 17;
130 s ^= s << 5; s >>>= 0;
131 return (s / 4294967296) * 2 - 1;
132 };
133 for (let k = 0; k < 2 * nlm; k++) q[k] = rnd();
134 for (let l = 0; l <= cfg.lmax; l++) q[2 * lmIndex(cfg.lmax, l, 0) + 1] = 0; // m=0 real
135 return q;
136}