Differentiate the phi flux in grid space
13 changed files+257−58
README.mdmodified+17−12View file
@@ -11,8 +11,9 @@ added is a *surface*.
1111
1212 The geometry is in the operator: the models evaluate the surface
1313 Laplace–Beltrami operator `lap_g` inside the implicit solve, in a **flux form
14-that costs 6 spherical-harmonic transforms per species per iteration** where
15-the textbook Cartesian-gradient form needs 12. See
14+that costs 5 spherical-harmonic transforms per species per iteration** (plus
15+one Legendre-free FFT derivative) where the textbook Cartesian-gradient form
16+needs 12. See
1617 [The geometry in the operator](#the-geometry-in-the-operator) and
1718 [docs/reduced-transforms.md](docs/reduced-transforms.md).
1819
@@ -115,11 +116,13 @@ for k = 1:niter
115116 % the sphere, one batched dispatch
116117 Pu = p1 .* Ftu + p2 .* Fpu; % the two fluxes, also smooth: the precomputed
117118 Qu = p2 .* Ftu + q2 .* Fpu; % weights carry every 1/sin(theta) there is
118- [PAu, QAu] = analys(Pu, Qu);
119+ PAu = analys(Pu);
119120 Pcu = PAu .* filt;
120- Qcu = QAu .* filt;
121- scu = dthetac(Pcu) + dphic(Qcu); % divergence, in coefficient space
122- lapu = r .* synth(scu); % = lap_g(u) on the grid
121+ scu = dthetac(Pcu); % theta part of the divergence, coefficients
122+ Lu = synth(scu); % sin(theta) * dtheta(P) on the grid
123+ dQu = dphig(Qu); % d/dphi is diagonal in the Fourier index:
124+ % two FFT stages, no Legendre work at all
125+ lapu = r .* (Lu + dQu); % = lap_g(u) on the grid
123126 dLu = (analys(lapu) + lamJ .* Un) .* filt; % dlap, projected onto the band
124127 Un = (Bu + (dt * D1) * dLu) ./ (1 + (dt * D1) * lamJ);
125128 end
@@ -169,9 +172,11 @@ Two formulations ship:
169172 ([`src/geom/metric.ts`](src/geom/metric.ts)), chosen so that **every field
170173 that gets analysed is a smooth function on the sphere** — the property
171174 that makes spherical-harmonic analysis meaningful, and the entire
172- difficulty near the poles. Cost: **6 transforms** per species per
173- iteration (3 syntheses + 3 analyses; `dthetac`/`dphic` are O(nlm)
174- coefficient shuffles, not transforms). The derivation, the smoothness
175+ difficulty near the poles. Cost: **5 Legendre transforms** per species
176+ per iteration (3 syntheses + 2 analyses; the phi flux never needs the
177+ Legendre basis — `dphig` differentiates it on the grid with two FFT
178+ stages, masking m past the top-degree filter — and `dthetac`/`dphic`
179+ are O(nlm) coefficient shuffles). The derivation, the smoothness
175180 argument and the fp32 error analysis are in
176181 [docs/reduced-transforms.md](docs/reduced-transforms.md).
177182 2. **The Cartesian-gradient form** (Algorithm 4 of `docs/algos.pdf`), kept as
@@ -424,9 +429,9 @@ Laplace-Beltrami scheme
424429 rests on;
425430 - on a non-axisymmetric surface, the flux tails match the Cartesian gradient
426431 component's, the doc's §7.1 criterion;
427-- the compiled op sequences add **6 transforms per species per iteration
428- against Algorithm 4's 12**, and a real simulation driven by each stays
429- within fp32 accumulation of the other.
432+- the compiled op sequences add **5 Legendre transforms per species per
433+ iteration against Algorithm 4's 12**, and a real simulation driven by
434+ each stays within fp32 accumulation of the other.
430435
431436 [`test/modelChecks.ts`](test/modelChecks.ts) compiles every model the app offers
432437 and asserts **how many kernels it compiles to**, split into the base step and
docs/reduced-transforms.mdmodified+13−0View file
@@ -128,6 +128,19 @@ Input $\{u^m_\ell\}$; output $\{(\Delta_\Gamma u)^m_\ell\}$.
128128 - The **only** division by $\sin\theta$ anywhere is folded into $p_1, p_2, q_2, r$ at precompute
129129 time. The per-matvec path contains none.
130130
131+### Implemented variation (2026-08-05): the $\varphi$-flux never needs the Legendre basis
132+
133+Step 4's analysis of $\tilde{Q}$ exists only so step 5 can apply $\partial_\varphi$ — but
134+$\partial_\varphi$ is diagonal in the Fourier index, so the implementation differentiates
135+$\tilde{Q}$ on the grid instead: FFT each latitude row, multiply mode $m$ by $im$ (zeroing
136+$m \ge L-2$ to mirror the top-degree filter; the Fourier analysis stage truncates $m > m_{\max}$
137+for free), inverse FFT. **5 Legendre transforms + one Legendre-free FFT derivative**, versus 6.
138+The caveat is that the grid route skips $\tilde Q$'s band projection in $\ell$; measured, this
139+does not bite — the band-edge spectra are identical to the 6-transform route's (the $m$ mask and
140+the final analysis's projection contain it), the Algorithm-4 A/B agreement is unchanged
141+($3.6\times10^{-4}$ after 20 steps at $L=63$), and the step gets ~8% faster at $L=255$
142+(~2% at $L=127$, where transform batching had already amortized most of what this removes).
143+
131144 ---
132145
133146 ## 5. Numerical trade-off
models/allencahn.mmodified+5−4View file
@@ -29,11 +29,12 @@ function [Un, u] = step(U, lam, filt, gx, gy, gz, p1, p2, q2, r, jhat, eps2, dt,
2929 [Ftu, Fpu] = synth(vtu, vpu);
3030 Pu = p1 .* Ftu + p2 .* Fpu;
3131 Qu = p2 .* Ftu + q2 .* Fpu;
32- [PAu, QAu] = analys(Pu, Qu);
32+ PAu = analys(Pu);
3333 Pcu = PAu .* filt;
34- Qcu = QAu .* filt;
35- scu = dthetac(Pcu) + dphic(Qcu);
36- lapu = r .* synth(scu);
34+ scu = dthetac(Pcu);
35+ Lu = synth(scu);
36+ dQu = dphig(Qu);
37+ lapu = r .* (Lu + dQu);
3738 dLu = (analys(lapu) + lamJ .* Un) .* filt;
3839
3940 Un = (Bu + (dt * eps2) * dLu) ./ (1 + (dt * eps2) * lamJ);
models/brusselator.mmodified+7−7View file
@@ -41,16 +41,16 @@ function [Un, Vn, u, v] = step(U, V, lam, filt, gx, gy, gz, p1, p2, q2, r, jhat,
4141 Qu = p2 .* Ftu + q2 .* Fpu;
4242 Pv = p1 .* Ftv + p2 .* Fpv;
4343 Qv = p2 .* Ftv + q2 .* Fpv;
44- [PAu, QAu, PAv, QAv] = analys(Pu, Qu, Pv, Qv);
44+ [PAu, PAv] = analys(Pu, Pv);
4545 Pcu = PAu .* filt;
46- Qcu = QAu .* filt;
4746 Pcv = PAv .* filt;
48- Qcv = QAv .* filt;
49- scu = dthetac(Pcu) + dphic(Qcu);
50- scv = dthetac(Pcv) + dphic(Qcv);
47+ scu = dthetac(Pcu);
48+ scv = dthetac(Pcv);
5149 [Lu, Lv] = synth(scu, scv);
52- lapu = r .* Lu;
53- lapv = r .* Lv;
50+ dQu = dphig(Qu);
51+ dQv = dphig(Qv);
52+ lapu = r .* (Lu + dQu);
53+ lapv = r .* (Lv + dQv);
5454 [LAu, LAv] = analys(lapu, lapv);
5555 dLu = (LAu + lamJ .* Un) .* filt;
5656 dLv = (LAv + lamJ .* Vn) .* filt;
models/schnakenberg.mmodified+16−13View file
@@ -58,10 +58,13 @@ function [Un, Vn, u, v] = step(U, V, lam, filt, gx, gy, gz, p1, p2, q2, r, jhat,
5858 % sin(theta)*dtheta(u) and dphi(u) -- both smooth on the sphere,
5959 % synthesized straight from the dthetac/dphic coefficient shuffles --
6060 % are combined pointwise through the precomputed weights p1,p2,q2 into
61- % two fluxes P,Q, also smooth. Their coefficients are then pushed
62- % through the *same* shuffles again and summed before the one synthesis
63- % of the divergence, which r scales into lap_g(u). The only division by
64- % sin(theta) anywhere is folded into p1,p2,q2,r at precompute time.
61+ % two fluxes P,Q, also smooth. The theta flux P goes back to
62+ % coefficients, through the same shuffle again, and is synthesized as
63+ % sin(theta)*dtheta(P); the phi flux Q never leaves the grid -- d/dphi
64+ % is diagonal in the Fourier index, so dphig differentiates it with two
65+ % FFT stages and no Legendre work (masking m past filt's reach). Their
66+ % sum, scaled by r, is lap_g(u). The only division by sin(theta)
67+ % anywhere is folded into p1,p2,q2,r at precompute time.
6568 % lamJ.*Un adds back the preconditioner's -lap_s(Un)/jhat, since lam
6669 % holds +l(l+1). filt zeroes the top two degrees, where the derivative
6770 % recurrences cannot exactly represent a derivative -- and the correction
@@ -72,8 +75,8 @@ function [Un, Vn, u, v] = step(U, V, lam, filt, gx, gy, gz, p1, p2, q2, r, jhat,
7275 % different rates -- a spurious Turing band at the band edge.
7376 %
7477 % The two species share each grouped call: the four gradient
75- % syntheses, the four flux analyses, the two divergence syntheses and
76- % the two final analyses each run as one batched dispatch.
78+ % syntheses, the two theta-flux analyses, the two divergence syntheses
79+ % and the two final analyses each run as one batched dispatch.
7780 Fu = Un .* filt;
7881 Fv = Vn .* filt;
7982 vtu = dthetac(Fu);
@@ -85,16 +88,16 @@ function [Un, Vn, u, v] = step(U, V, lam, filt, gx, gy, gz, p1, p2, q2, r, jhat,
8588 Qu = p2 .* Ftu + q2 .* Fpu;
8689 Pv = p1 .* Ftv + p2 .* Fpv;
8790 Qv = p2 .* Ftv + q2 .* Fpv;
88- [PAu, QAu, PAv, QAv] = analys(Pu, Qu, Pv, Qv);
91+ [PAu, PAv] = analys(Pu, Pv);
8992 Pcu = PAu .* filt;
90- Qcu = QAu .* filt;
9193 Pcv = PAv .* filt;
92- Qcv = QAv .* filt;
93- scu = dthetac(Pcu) + dphic(Qcu);
94- scv = dthetac(Pcv) + dphic(Qcv);
94+ scu = dthetac(Pcu);
95+ scv = dthetac(Pcv);
9596 [Lu, Lv] = synth(scu, scv);
96- lapu = r .* Lu;
97- lapv = r .* Lv;
97+ dQu = dphig(Qu);
98+ dQv = dphig(Qv);
99+ lapu = r .* (Lu + dQu);
100+ lapv = r .* (Lv + dQv);
98101 [LAu, LAv] = analys(lapu, lapv);
99102 dLu = (LAu + lamJ .* Un) .* filt;
100103 dLv = (LAv + lamJ .* Vn) .* filt;
src/mgpu/externals.tsmodified+9−1View file
@@ -145,10 +145,18 @@ export function externalOpFiles(g: GridSizes): { name: string; source: string }[
145145 name: 'dphic.mtoc2.js',
146146 source: transformSource('dphic', 2, g.nlm, 2, g.nlm),
147147 },
148+ {
149+ // Grid-space phi-derivative: two Fourier stages and an i*m multiply,
150+ // no Legendre work (d/dphi is diagonal in the Fourier index). What
151+ // lets the flux-form divergence skip the Q-flux's spherical-harmonic
152+ // analysis.
153+ name: 'dphig.mtoc2.js',
154+ source: transformSource('dphig', g.npts, 1, g.npts, 1),
155+ },
148156 ];
149157 }
150158
151159 /** Names the WGSL backend must implement as GPU encodes rather than kernels. */
152160 export const EXTERNAL_OPS = new Set([
153- 'synth', 'analys', 'dtheta', 'dphi', 'dthetac', 'dphic',
161+ 'synth', 'analys', 'dtheta', 'dphi', 'dthetac', 'dphic', 'dphig',
154162 ]);
src/mgpu/plan.tsmodified+15−1View file
@@ -17,7 +17,7 @@ import type {
1717 MultiAssignCall,
1818 } from 'numbl-src/numbl-core/jit/lowering/ir.ts';
1919 import type { NumericType, Type } from 'numbl-src/numbl-core/jit/lowering/types.ts';
20-import { ShtPlan, type ShtBinding, type ShtBatchBinding } from '../sht/sht.ts';
20+import { ShtPlan, type ShtBinding, type ShtBatchBinding, type ShtDphigBinding } from '../sht/sht.ts';
2121 import { DerivPlan, type DerivBinding } from '../sht/deriv.ts';
2222 import type { CompiledFunction } from './compile.ts';
2323 import { EXTERNAL_OPS } from './externals.ts';
@@ -129,6 +129,7 @@ type Op =
129129 | { kind: 'synth-batch' | 'analys-batch'; binding: ShtBatchBinding; labels: string[] }
130130 | { kind: 'dtheta' | 'dphi'; binding: DerivBinding; label: string }
131131 | { kind: 'dthetac' | 'dphic'; bindGroup: GPUBindGroup; label: string }
132+ | { kind: 'dphig'; binding: ShtDphigBinding; label: string }
132133 | { kind: 'copy'; from: GPUBuffer; to: GPUBuffer; bytes: number; label: string };
133134
134135 /**
@@ -464,6 +465,16 @@ export class ModelPlan {
464465 );
465466 }
466467 const label = `${stmt.name} = ${ext.name}(${ext.argName})`;
468+ if (ext.name === 'dphig') {
469+ // Grid -> grid, staged through the plan's fm scratch; safe even
470+ // in place, so no aliasing guard is needed.
471+ planned.push({
472+ kind: 'dphig',
473+ binding: sht.createDphigBinding(argSlot.buffer, dest.buffer),
474+ label,
475+ });
476+ return;
477+ }
467478 if (ext.name === 'synth' || ext.name === 'analys') {
468479 // Left unbound until materializeTransforms has grouped adjacent
469480 // independent transforms into batched dispatches.
@@ -779,6 +790,9 @@ export class ModelPlan {
779790 case 'dphic':
780791 this.#deriv!.encodeDphicInto(inPass(), op.bindGroup);
781792 break;
793+ case 'dphig':
794+ this.#sht.encodeDphigInto(inPass(), op.binding);
795+ break;
782796 case 'synth-batch':
783797 this.#sht.encodeSynthBatchInto(inPass(), op.binding);
784798 break;
src/sht/sht.tsmodified+87−1View file
@@ -15,6 +15,7 @@ import {
1515 legSynthBatchWGSL,
1616 legAnalysBatchWGSL,
1717 } from './wgsl/leg.ts';
18+import { fmDphiWGSL } from './wgsl/deriv.ts';
1819 import {
1920 fftSynthWGSL,
2021 fftAnalysWGSL,
@@ -33,6 +34,13 @@ export interface ShtBinding {
3334 readonly bgFour: GPUBindGroup;
3435 }
3536
37+/** The three bind groups of one grid-space phi-derivative (see dphig). */
38+export interface ShtDphigBinding {
39+ readonly bgFourAnalys: GPUBindGroup;
40+ readonly bgMul: GPUBindGroup;
41+ readonly bgFourSynth: GPUBindGroup;
42+}
43+
3644 /**
3745 * One batched transform: K fields through a single Legendre dispatch (the
3846 * recurrence walked once, K accumulator lanes) plus K per-field Fourier
@@ -175,6 +183,8 @@ export class ShtPlan {
175183 private pipeLegAnalys!: GPUComputePipeline;
176184 private pipeFourSynth!: GPUComputePipeline;
177185 private pipeFourAnalys!: GPUComputePipeline;
186+ /** Fourier-space i*m multiply, the middle of dphig. */
187+ private pipeFmDphi!: GPUComputePipeline;
178188 /** Batched Legendre pipelines by lane count (even sizes up to batchK). */
179189 private pipeLegSynthB = new Map<number, GPUComputePipeline>();
180190 private pipeLegAnalysB = new Map<number, GPUComputePipeline>();
@@ -289,7 +299,7 @@ export class ShtPlan {
289299 this.fourierMode === 'fft' && nphi % 2 === 0 && tuning('SHT_REAL_FFT') !== false;
290300 const fftS = realFft ? fftSynthRealWGSL : fftSynthWGSL;
291301 const fftA = realFft ? fftAnalysRealWGSL : fftAnalysWGSL;
292- const [pLegS, pLegA, pFourS, pFourA] = await Promise.all([
302+ const [pLegS, pLegA, pFourS, pFourA, pFmDphi] = await Promise.all([
293303 makePipeline(dev, legSynthWGSL(legP), 'leg_synth'),
294304 makePipeline(dev, legAnalysWGSL(legP), 'leg_analys'),
295305 makePipeline(
@@ -302,11 +312,19 @@ export class ShtPlan {
302312 this.fourierMode === 'fft' ? fftA(fourP) : dftAnalysWGSL(fourP),
303313 this.fourierMode === 'fft' ? 'fft_analys' : 'dft_analys',
304314 ),
315+ // Mirrors filterMask: content at l >= lmax-2 is filtered on the
316+ // l-space route, so the m-space route keeps m <= lmax-3.
317+ makePipeline(
318+ dev,
319+ fmDphiWGSL({ mmax, nlat, nphi, mcut: lmax - 3 }),
320+ 'fm_dphi',
321+ ),
305322 ]);
306323 this.pipeLegSynth = pLegS;
307324 this.pipeLegAnalys = pLegA;
308325 this.pipeFourSynth = pFourS;
309326 this.pipeFourAnalys = pFourA;
327+ this.pipeFmDphi = pFmDphi;
310328
311329 // --- batched Legendre pipelines ---
312330 // The widest even K <= 4 whose bind group (3 tables + K fields + the fm
@@ -472,6 +490,74 @@ export class ShtPlan {
472490 };
473491 }
474492
493+ /**
494+ * Bind groups for one grid-space phi-derivative, dphig: Fourier analysis
495+ * of each latitude row into fm (which truncates to m <= mmax for free),
496+ * the i*m/NPHI multiply (zeroing m past the top-degree filt's reach), and
497+ * Fourier synthesis back to the grid. No Legendre stage anywhere — this
498+ * is what lets the flux-form divergence drop the Q-flux's spherical-
499+ * harmonic analysis (docs/reduced-transforms.md Sec 5b's companion trick
500+ * in Sec 6-of-changes): d/dphi is diagonal in the Fourier index. Uses
501+ * fmBuf as scratch, sequentially like every transform in a pass.
502+ */
503+ createDphigBinding(spatIn: GPUBuffer, spatOut: GPUBuffer): ShtDphigBinding {
504+ return {
505+ bgFourAnalys: this.device.createBindGroup({
506+ layout: this.pipeFourAnalys.getBindGroupLayout(0),
507+ entries: bgEntries([spatIn, this.fmBuf, this.bufTrig]),
508+ }),
509+ bgMul: this.device.createBindGroup({
510+ layout: this.pipeFmDphi.getBindGroupLayout(0),
511+ entries: bgEntries([this.fmBuf]),
512+ }),
513+ bgFourSynth: this.device.createBindGroup({
514+ layout: this.pipeFourSynth.getBindGroupLayout(0),
515+ entries: bgEntries([this.fmBuf, spatOut, this.bufTrig]),
516+ }),
517+ };
518+ }
519+
520+ /** Record dphig into an existing compute pass: two Fourier stages and a
521+ * pointwise multiply — no Legendre work. */
522+ encodeDphigInto(pass: GPUComputePassEncoder, b: ShtDphigBinding): void {
523+ const { mmax, nlat, nphi } = this.cfg;
524+ pass.setPipeline(this.pipeFourAnalys);
525+ pass.setBindGroup(0, b.bgFourAnalys);
526+ if (this.fourierMode === 'fft') {
527+ pass.dispatchWorkgroups(nlat);
528+ } else {
529+ pass.dispatchWorkgroups(Math.ceil((mmax + 1) / 64), nlat);
530+ }
531+ pass.setPipeline(this.pipeFmDphi);
532+ pass.setBindGroup(0, b.bgMul);
533+ pass.dispatchWorkgroups(Math.ceil(((mmax + 1) * nlat) / 64));
534+ pass.setPipeline(this.pipeFourSynth);
535+ pass.setBindGroup(0, b.bgFourSynth);
536+ if (this.fourierMode === 'fft') {
537+ pass.dispatchWorkgroups(nlat);
538+ } else {
539+ pass.dispatchWorkgroups(Math.ceil(nphi / 64), nlat);
540+ }
541+ }
542+
543+ /** CPU convenience: grid field -> d/dphi of its trig interpolant, for tests. */
544+ async dphig(spat: Float32Array): Promise<Float32Array> {
545+ const { nlat, nphi } = this.cfg;
546+ if (spat.length !== nlat * nphi) throw new Error(`spat must have length ${nlat * nphi}`);
547+ this.device.queue.writeBuffer(this.spatBuf, 0, spat as Float32Array<ArrayBuffer>);
548+ const binding = this.createDphigBinding(this.spatBuf, this.spatBuf);
549+ const enc = this.device.createCommandEncoder({ label: 'sht-dphig' });
550+ const pass = enc.beginComputePass({ label: 'sht-dphig' });
551+ this.encodeDphigInto(pass, binding);
552+ pass.end();
553+ enc.copyBufferToBuffer(this.spatBuf, 0, this.stageSpat, 0, 4 * nlat * nphi);
554+ this.device.queue.submit([enc.finish()]);
555+ await this.stageSpat.mapAsync(GPUMapMode.READ);
556+ const out = new Float32Array(this.stageSpat.getMappedRange().slice(0));
557+ this.stageSpat.unmap();
558+ return out;
559+ }
560+
475561 /** Record a batched synthesis: one Legendre dispatch, K Fourier dispatches. */
476562 encodeSynthBatchInto(pass: GPUComputePassEncoder, b: ShtBatchBinding): void {
477563 const { mmax, nlat, nphi } = this.cfg;
src/sht/wgsl/deriv.tsmodified+41−0View file
@@ -99,3 +99,44 @@ fn divide_sin_theta(@builtin(global_invocation_id) gid: vec3u) {
9999 }
100100 `;
101101 }
102+
103+export interface FmDphiParams {
104+ mmax: number;
105+ nlat: number;
106+ nphi: number;
107+ /** Highest m kept; modes above are zeroed (mirrors the l-space filt). */
108+ mcut: number;
109+}
110+
111+/**
112+ * The Fourier-space middle of the grid-space phi-derivative `dphig`:
113+ * fm holds the unnormalized DFT modes of each latitude row (what the
114+ * Fourier analysis stage produces), so d/dphi is fm[m] *= i*m/NPHI --
115+ * the 1/NPHI undoes the unnormalized analysis+synthesis round trip.
116+ * Modes above MCUT are zeroed: the Fourier analysis stage already
117+ * truncated m > mmax for free, and MCUT additionally mirrors the
118+ * top-degree filt so the differentiated field carries no content the
119+ * l-space route would not have kept.
120+ */
121+export function fmDphiWGSL(p: FmDphiParams): string {
122+ const count = (p.mmax + 1) * p.nlat;
123+ return /* wgsl */ `
124+const NLAT: u32 = ${p.nlat}u;
125+const COUNT: u32 = ${count}u;
126+const MCUT: u32 = ${Math.max(0, p.mcut)}u;
127+const INV_NPHI: f32 = ${1 / p.nphi};
128+
129+@group(0) @binding(0) var<storage, read_write> fm: array<vec2f>;
130+
131+@compute @workgroup_size(${WG})
132+fn fm_dphi(@builtin(global_invocation_id) gid: vec3u) {
133+ let i = gid.x;
134+ if (i >= COUNT) { return; }
135+ let m = i / NLAT;
136+ var k: f32 = 0.0;
137+ if (m <= MCUT) { k = f32(m) * INV_NPHI; }
138+ let c = fm[i];
139+ fm[i] = vec2f(-k * c.y, k * c.x);
140+}
141+`;
142+}
test/fluxChecks.tsmodified+6−4View file
@@ -368,11 +368,13 @@ export async function fluxChecks(
368368 xformsPerIter.push(counts[1] - counts[0]);
369369 }
370370
371- // The headline number, from the compiled op sequences: 6 transforms per
372- // species per iteration against Algorithm 4's 12 (2 species here).
371+ // The headline number, from the compiled op sequences: 5 Legendre
372+ // transforms per species per iteration against Algorithm 4's 12
373+ // (2 species here). The phi flux's derivative runs as dphig -- two
374+ // Fourier stages, no Legendre work -- and is deliberately not counted.
373375 check(
374- 'flux: 6 transforms per species per iteration, versus 12',
375- xformsPerIter[0] === 12 && xformsPerIter[1] === 24,
376+ 'flux: 5 Legendre transforms per species per iteration, versus 12',
377+ xformsPerIter[0] === 10 && xformsPerIter[1] === 24,
376378 `flux form adds ${xformsPerIter[0]} transforms/iteration, ` +
377379 `Algorithm 4 adds ${xformsPerIter[1]}`,
378380 );
test/geometryChecks.tsmodified+6−6View file
@@ -311,13 +311,13 @@ export async function geometryChecks(
311311 );
312312 // Unrolling has to be exactly linear in the trip count: the body planned
313313 // once per iteration, no more and no less. Per species per iteration: 3
314- // synths + 3 analyses (the flux-form matvec's six real transforms,
315- // docs/reduced-transforms.md Sec 4) + 4 coefficient-space
316- // dthetac/dphic shuffles plus 9 generated kernels -- see
317- // test/modelChecks.ts's KERNELS_PER_ITERATION, which counts the kernels
318- // alone; this counts every op, transforms and shuffles included.
314+ // synths + 2 analyses (the flux-form matvec's five Legendre transforms,
315+ // docs/reduced-transforms.md Sec 4 with the dphig variation) + the
316+ // grid-space phi-derivative + 3 coefficient-space shuffles plus 7
317+ // generated kernels -- see test/modelChecks.ts's KERNELS_PER_ITERATION,
318+ // which counts the kernels alone; this counts every op.
319319 const perIteration = ops[1] - ops[0];
320- const want = 38;
320+ const want = 32;
321321 check(
322322 'loop: unrolling is exactly linear in the trip count',
323323 perIteration === want && ops[2] - ops[0] === 4 * perIteration,
test/modelChecks.tsmodified+10−9View file
@@ -47,9 +47,9 @@ const EXPECTED_KERNELS: Record<string, number> = {
4747 * as a live reference, with its original counts.
4848 */
4949 const KERNELS_PER_ITERATION: Record<string, number> = {
50- schnakenberg: 18,
51- brusselator: 18,
52- allencahn: 9,
50+ schnakenberg: 14,
51+ brusselator: 14,
52+ allencahn: 7,
5353 // 30 before the correction gained its band projection (.* filt on dLu):
5454 // that line fused into the state update in this model's expression shape,
5555 // and no longer does — one extra 2 x nlm kernel per species per iteration.
@@ -178,14 +178,15 @@ export async function modelChecks(
178178 }
179179 // Every batchable run at one solve iteration: the u/v syntheses and the
180180 // reaction analyses outside the loop (2 + 2), the four gradient
181- // syntheses, four flux analyses, two divergence syntheses and two final
182- // analyses inside it (4 + 4 + 2 + 2). Lane counts are batch-width
183- // invariant: a x4 run is one batch at K = 4 and two at K = 2, but the
184- // lanes annotated are the same 16 either way.
181+ // syntheses, two theta-flux analyses, two divergence syntheses and two
182+ // final analyses inside it (4 + 2 + 2 + 2; the phi flux goes through
183+ // dphig, which has no Legendre stage to batch). Lane counts are
184+ // batch-width invariant: a x4 run is one batch at K = 4 and two at
185+ // K = 2, but the lanes annotated are the same 14 either way.
185186 check(
186187 'batch: the compiled step batches every adjacent transform pair',
187- batchedLanes === 16,
188- `${batchedLanes} batched transform lanes (expected 16)`,
188+ batchedLanes === 14,
189+ `${batchedLanes} batched transform lanes (expected 14)`,
189190 );
190191 let worst = 0;
191192 for (let i = 0; i < states[0].length; i++) {
test/transformChecks.tsmodified+25−0View file
@@ -126,6 +126,31 @@ export async function transformChecks(
126126 deriv.destroy();
127127 }
128128
129+ // ---- grid-space phi-derivative (dphig) vs the f64 reference -------------
130+ // dphig differentiates in phi with two Fourier stages and an i*m multiply,
131+ // no Legendre work. On a band-limited field whose m >= lmax-2 modes are
132+ // zero (dphig masks those, mirroring filt), it must agree with the
133+ // coefficient-space route dphi = synth(i*m*coeffs) to fp32.
134+ {
135+ const q64 = new Float64Array(randomSpectrum(cfg, 4242));
136+ for (let m = Math.max(0, lmax - 2); m <= lmax; m++) {
137+ for (let l = m; l <= lmax; l++) {
138+ const i = 2 * (m * (lmax + 1) - (m * (m - 1)) / 2 + (l - m));
139+ q64[i] = 0;
140+ q64[i + 1] = 0;
141+ }
142+ }
143+ const grid = ref.synth(q64);
144+ const dPhiGpu = await plan.dphig(new Float32Array(grid));
145+ const dPhiCpu = ref.dphi(q64);
146+ const err = relL2(dPhiGpu, dPhiCpu);
147+ check(
148+ 'dphig: grid-space FFT phi-derivative vs f64 CPU reference',
149+ err < 1e-4,
150+ `rel L2 ${err.toExponential(2)}`,
151+ );
152+ }
153+
129154 // ---- batched transforms reproduce the scalar transforms ------------------
130155 // A batch walks the Legendre recurrence once for K fields with per-lane
131156 // arithmetic textually identical to the scalar kernel's, so each lane must