forked from github.com/concept-collection/turing-surface
2 * Grid and spectral layout definitions, following SHTNS conventions:
3 *
4 * - Spectral coefficients Q_lm are complex, stored for m >= 0 only (real
5 * fields), interleaved [re, im], with SHTNS "m-major" ordering:
6 * for m = 0..mmax: for l = m..lmax. Index of (l, m) is lm(l, m).
7 * - Spatial fields are real, phi-contiguous: spat[ilat * nphi + iphi],
8 * with ilat ordered by increasing colatitude theta (north to south)
9 * and iphi covering [0, 2*pi) uniformly.
10 * - Normalization: orthonormal spherical harmonics INCLUDING the
11 * Condon-Shortley phase (SHTNS default: sht_orthonormal).
12 * A real field is f = sum_{l,m>=0} Q_lm Y_lm + c.c.(m>0), i.e.
13 * Q_{l,-m} = (-1)^m conj(Q_lm) is implied. m=0 coefficients must
14 * have zero imaginary part.
15 */
17export interface ShtConfig {
18 lmax: number;
19 mmax: number;
20 nlat: number;
21 nphi: number;
22}
24export function nlmCalc(lmax: number, mmax: number): number {
25 // sum over m=0..mmax of (lmax - m + 1)
26 return (mmax + 1) * (lmax + 1) - (mmax * (mmax + 1)) / 2;
27}
29/** Index of coefficient (l, m) in the spectral array (SHTNS LM ordering). */
30export function lmIndex(lmax: number, l: number, m: number): number {
31 return m * (lmax + 1) - (m * (m - 1)) / 2 + (l - m);
32}
34export function validateConfig(cfg: ShtConfig): void {
35 const { lmax, mmax, nlat, nphi } = cfg;
36 if (!Number.isInteger(lmax) || lmax < 1) throw new Error(`lmax must be an integer >= 1 (got ${lmax})`);
37 if (!Number.isInteger(mmax) || mmax < 0 || mmax > lmax)
38 throw new Error(`mmax must be an integer in [0, lmax] (got ${mmax})`);
39 if (!Number.isInteger(nlat) || nlat <= lmax)
40 throw new Error(`nlat must be an integer > lmax for exact Gauss quadrature (got nlat=${nlat}, lmax=${lmax})`);
41 if (!Number.isInteger(nphi) || nphi < 2 * mmax + 1)
42 throw new Error(`nphi must be an integer >= 2*mmax+1 to avoid aliasing (got nphi=${nphi}, mmax=${mmax})`);
43}
45export function isPowerOfTwo(n: number): boolean {
46 return n > 0 && (n & (n - 1)) === 0;
47}
50 * Re-index coefficients from one (lmax, mmax) layout into another, dropping
51 * the modes the target cannot hold and zero-filling the ones the source does
52 * not have. Both truncation and padding are exact spectral operations —
53 * truncation is orthogonal projection onto the smaller band, padding changes
54 * nothing — so a surface carried as coefficients can be moved onto any plan.
55 */
56export function relayoutCoeffs(
57 q: Float32Array,
58 from: { lmax: number; mmax: number },
59 to: { lmax: number; mmax: number },
60): Float32Array {
61 if (from.lmax === to.lmax && from.mmax === to.mmax) return q;
62 const out = new Float32Array(2 * nlmCalc(to.lmax, to.mmax));
63 const mTop = Math.min(from.mmax, to.mmax);
64 const lTop = Math.min(from.lmax, to.lmax);
65 for (let m = 0; m <= mTop; m++) {
66 for (let l = m; l <= lTop; l++) {
67 const src = 2 * lmIndex(from.lmax, l, m);
68 const dst = 2 * lmIndex(to.lmax, l, m);
69 out[dst] = q[src];
70 out[dst + 1] = q[src + 1];
71 }
72 }
73 return out;
74}
0af3386Reaction-diffusion on spherical-harmonic surfacesJeremy Magland 76/** Grid sizes for a given lmax, dealiased for a reaction of polynomial degree
77 * `pdeg` (the rule from websph's reference implementation):
78 * nlat >= ((pdeg+1)*lmax+1)/2, nphi >= (pdeg+1)*lmax+1. nphi is rounded up to
79 * a power of two to keep the GPU FFT path. */
80export function gridForLmax(lmax: number, pdeg: number): { nlat: number; nphi: number } {
81 const minLat = Math.max(lmax + 1, ((pdeg + 1) * lmax + 1) / 2);
82 const nlat = 2 * Math.ceil(minLat / 2);
83 let nphi = 1;
84 while (nphi < (pdeg + 1) * lmax + 1) nphi *= 2;
85 return { nlat, nphi };
86}