2dedc35turing-sphere: reaction-diffusion on the sphere, spectral solver on WebGPUJeremy Magland 1/**
2 * WebGPU spherical harmonic transform plan (scalar transforms, fp32).
3 *
4 * Mirrors the structure of the SHTNS CUDA backend (sht_gpu.cu):
5 * host-side f64 precomputation of grid + recurrence coefficients, shader
6 * source generated with sizes baked in (SHTNS uses NVRTC; WGSL is always
7 * runtime-compiled), then per-transform: Legendre stage + Fourier stage.
8 */
9import { gaussNodesWeights } from './gauss.ts';
10import { legendreCoeffs } from './coeffs.ts';
11import { nlmCalc, validateConfig, isPowerOfTwo, type ShtConfig } from './layout.ts';
12import { legSynthWGSL, legAnalysWGSL } from './wgsl/leg.ts';
13import {
14 fftSynthWGSL,
15 fftAnalysWGSL,
16 dftSynthWGSL,
17 dftAnalysWGSL,
18 fftThreads,
19} from './wgsl/fourier.ts';
21export type FourierMode = 'auto' | 'fft' | 'dft';
61e12f1Write the solver in MATLAB and compile it to WebGPUJeremy Magland 23/** The two bind groups (Legendre stage, Fourier stage) of one transform. */
24export interface ShtBinding {
25 readonly bgLeg: GPUBindGroup;
26 readonly bgFour: GPUBindGroup;
27}
29const bgEntries = (bufs: GPUBuffer[]) =>
30 bufs.map((buffer, binding) => ({ binding, resource: { buffer } }));
2dedc35turing-sphere: reaction-diffusion on the sphere, spectral solver on WebGPUJeremy Magland 32export interface ShtOptions {
33 /** Fourier stage implementation. 'auto' picks fft when nphi is a power of two that fits in workgroup memory. */
34 fourier?: FourierMode;
35}
37const WG_SYNTH = 64;
39/**
40 * Tuning knob, for A/B-ing a change without editing code. Reads globalThis
41 * first (set it before creating a plan, as scripts/_ab.ts does), then the
42 * environment, so `SHT_SUBGROUPS=0 npm run bench:sht` works too. `process` is
43 * absent in the browser, where only the globalThis form applies.
44 */
45function tuning(name: string): unknown {
46 const g = (globalThis as Record<string, unknown>)[name];
47 if (g !== undefined) return g;
48 const env = (globalThis as { process?: { env?: Record<string, string | undefined> } }).process?.env?.[name];
49 if (env === undefined || env === '') return undefined;
50 if (env === '1' || env === 'true') return true;
51 if (env === '0' || env === 'false') return false;
52 const n = Number(env);
53 return Number.isFinite(n) ? n : env;
54}
56/**
57 * Workgroup size for the analysis Legendre reduction. The right answer differs
58 * between the two reduction strategies, so it is chosen per strategy.
59 *
60 * Measured on an RTX PRO 6000 Blackwell (analysis, us). Shared-memory tree,
61 * where each doubling of wgAnalys costs another barrier per l-pair:
62 *
63 * wgAnalys: 16 32 64 128 256
64 * nlat=128 43.2 40.6 42.6 46.7 52.9 -> 32
65 * nlat=256 104.0 85.6 87.8 89.9 100.0 -> 32
66 * nlat=512 408.8 216.9 184.4 190.8 204.8 -> 64
67 *
68 * i.e. max(32, nlat/8). A flat 32 would be worse than the old default of 256 at
69 * nlat=512, so it cannot be fitted on one grid. With subgroupAdd the barrier
70 * count stops growing with wgAnalys and the picture inverts: threads in flight,
71 * (mmax+1) * wgAnalys, becomes binding, since analysis dispatches only mmax+1
72 * workgroups. 128 then wins at every grid (round trip, us):
73 *
74 * 128x256 37.8 (vs 38.3), 256x512 66.6 (vs 74.7), 512x1024 131.9 (vs 156.6)
75 */
76function defaultWgAnalys(nlat: number, limit: number, subgroups: boolean): number {
77 if (subgroups) {
78 // capped at nlat so small grids do not launch threads with no latitude to own
79 let cap = 1;
80 while (cap < nlat) cap *= 2;
81 return Math.min(128, limit, cap);
82 }
83 const target = Math.max(32, nlat / 8);
84 let wg = 1;
85 while (wg < target) wg *= 2; // the tree reduction halves, so a power of two
86 return Math.min(wg, limit);
87}
89async function makePipeline(
90 device: GPUDevice,
91 code: string,
92 entryPoint: string,
93): Promise<GPUComputePipeline> {
94 device.pushErrorScope('validation');
95 const module = device.createShaderModule({ code, label: entryPoint });
96 const info = await module.getCompilationInfo();
97 const errors = info.messages.filter((m) => m.type === 'error');
98 if (errors.length) {
99 throw new Error(
100 `WGSL compile error in ${entryPoint}:\n` +
101 errors.map((e) => ` ${e.lineNum}:${e.linePos} ${e.message}`).join('\n'),
102 );
103 }
104 const pipeline = await device.createComputePipelineAsync({
105 layout: 'auto',
106 compute: { module, entryPoint },
107 label: entryPoint,
108 });
109 const err = await device.popErrorScope();
110 if (err) throw new Error(`pipeline ${entryPoint}: ${err.message}`);
111 return pipeline;
112}
114export class ShtPlan {
115 readonly cfg: ShtConfig;
116 readonly nlm: number;
117 readonly fourierMode: 'fft' | 'dft';
eb4a9e7Fold north/south latitude pairs onto one Legendre recurrencedanfortunato 118 /** Latitudes leg_synth walks: nlat/2 when parity folding. */
119 readonly legLat: number = 0;
2dedc35turing-sphere: reaction-diffusion on the sphere, spectral solver on WebGPUJeremy Magland 120 /** Colatitudes theta_i (f64, increasing: north to south). */
121 readonly theta: Float64Array;
122 readonly cosTheta: Float64Array;
123 readonly gaussWeights: Float64Array;
125 private device: GPUDevice;
126 private bufAb!: GPUBuffer;
127 private bufAmm!: GPUBuffer;
128 private bufCtstw!: GPUBuffer;
129 private bufTrig!: GPUBuffer;
130 /** Spectral input (synthesis) — write with queue.writeBuffer or use synth(). */
131 readonly qlmIn!: GPUBuffer;
132 /** Spectral output (analysis). */
133 readonly qlmOut!: GPUBuffer;
30ed90fAdd scripts/diagnose-sht.ts: which stage of the transform is wrong?Jeremy Magland 134 /** Fourier-space intermediate [(m)*nlat + ilat], complex f32. COPY_SRC so the
135 * stage boundary is observable: a transform is Legendre-then-Fourier, and
136 * scripts/diagnose-sht.ts tells the two apart by reading this. */
2dedc35turing-sphere: reaction-diffusion on the sphere, spectral solver on WebGPUJeremy Magland 137 readonly fmBuf!: GPUBuffer;
138 /** Spatial field [ilat*nphi + iphi], f32. */
139 readonly spatBuf!: GPUBuffer;
140 private stageSpat!: GPUBuffer;
141 private stageQ!: GPUBuffer;
143 private pipeLegSynth!: GPUComputePipeline;
144 private pipeLegAnalys!: GPUComputePipeline;
145 private pipeFourSynth!: GPUComputePipeline;
146 private pipeFourAnalys!: GPUComputePipeline;
147 private bgLegSynth!: GPUBindGroup;
148 private bgLegAnalys!: GPUBindGroup;
149 private bgFourSynth!: GPUBindGroup;
150 private bgFourAnalys!: GPUBindGroup;
152 private constructor(device: GPUDevice, cfg: ShtConfig, fourierMode: 'fft' | 'dft') {
153 this.device = device;
154 this.cfg = cfg;
155 this.nlm = nlmCalc(cfg.lmax, cfg.mmax);
156 this.fourierMode = fourierMode;
157 const { x, w } = gaussNodesWeights(cfg.nlat);
158 this.cosTheta = x;
159 this.gaussWeights = w;
160 this.theta = new Float64Array(cfg.nlat);
161 for (let i = 0; i < cfg.nlat; i++) this.theta[i] = Math.acos(x[i]);
162 }
164 static async create(device: GPUDevice, cfg: ShtConfig, opts: ShtOptions = {}): Promise<ShtPlan> {
165 validateConfig(cfg);
166 const want = opts.fourier ?? 'auto';
167 const fftFits =
168 isPowerOfTwo(cfg.nphi) &&
169 16 * cfg.nphi <= device.limits.maxComputeWorkgroupStorageSize &&
170 fftThreads(cfg.nphi) <= device.limits.maxComputeInvocationsPerWorkgroup;
171 if (want === 'fft' && !fftFits) {
172 throw new Error(
173 `fourier:'fft' requires power-of-two nphi with 16*nphi <= maxComputeWorkgroupStorageSize ` +
174 `(nphi=${cfg.nphi}, limit=${device.limits.maxComputeWorkgroupStorageSize})`,
175 );
176 }
177 const mode: 'fft' | 'dft' = want === 'dft' ? 'dft' : fftFits ? 'fft' : 'dft';
178 const plan = new ShtPlan(device, cfg, mode);
179 await plan.init();
180 return plan;
181 }
183 private async init(): Promise<void> {
184 const { lmax, mmax, nlat, nphi } = this.cfg;
185 const dev = this.device;
186 const self = this as {
187 -readonly [k in keyof ShtPlan]: ShtPlan[k];
188 };
190 // --- host precomputation (f64), then downcast to f32 for upload ---
191 const { amm, ab } = legendreCoeffs(lmax, mmax);
192 const ctstw = new Float32Array(3 * nlat);
193 for (let i = 0; i < nlat; i++) {
194 ctstw[i] = this.cosTheta[i];
195 ctstw[nlat + i] = Math.sqrt(1 - this.cosTheta[i] * this.cosTheta[i]);
196 ctstw[2 * nlat + i] = this.gaussWeights[i] * ((2 * Math.PI) / nphi);
197 }
198 // twiddle/phase table in f64 (device sin/cos is too inaccurate: ~2^-11 under Vulkan)
199 const trig = new Float32Array(2 * nphi);
200 for (let k = 0; k < nphi; k++) {
201 trig[2 * k] = Math.cos((2 * Math.PI * k) / nphi);
202 trig[2 * k + 1] = Math.sin((2 * Math.PI * k) / nphi);
203 }
205 const mkBuf = (label: string, size: number, usage: GPUBufferUsageFlags) =>
206 dev.createBuffer({ label, size, usage });
207 this.bufAb = mkBuf('sht-ab', 8 * this.nlm, GPUBufferUsage.STORAGE | GPUBufferUsage.COPY_DST);
208 this.bufAmm = mkBuf('sht-amm', 4 * (mmax + 1), GPUBufferUsage.STORAGE | GPUBufferUsage.COPY_DST);
209 this.bufCtstw = mkBuf('sht-ctstw', 4 * 3 * nlat, GPUBufferUsage.STORAGE | GPUBufferUsage.COPY_DST);
210 this.bufTrig = mkBuf('sht-trig', 8 * nphi, GPUBufferUsage.STORAGE | GPUBufferUsage.COPY_DST);
211 self.qlmIn = mkBuf('sht-qlm-in', 8 * this.nlm, GPUBufferUsage.STORAGE | GPUBufferUsage.COPY_DST | GPUBufferUsage.COPY_SRC);
212 self.qlmOut = mkBuf('sht-qlm-out', 8 * this.nlm, GPUBufferUsage.STORAGE | GPUBufferUsage.COPY_SRC);
30ed90fAdd scripts/diagnose-sht.ts: which stage of the transform is wrong?Jeremy Magland 213 self.fmBuf = mkBuf('sht-fm', 8 * (mmax + 1) * nlat, GPUBufferUsage.STORAGE | GPUBufferUsage.COPY_SRC);
2dedc35turing-sphere: reaction-diffusion on the sphere, spectral solver on WebGPUJeremy Magland 214 self.spatBuf = mkBuf('sht-spat', 4 * nlat * nphi, GPUBufferUsage.STORAGE | GPUBufferUsage.COPY_DST | GPUBufferUsage.COPY_SRC);
215 this.stageSpat = mkBuf('sht-stage-spat', 4 * nlat * nphi, GPUBufferUsage.MAP_READ | GPUBufferUsage.COPY_DST);
216 this.stageQ = mkBuf('sht-stage-q', 8 * this.nlm, GPUBufferUsage.MAP_READ | GPUBufferUsage.COPY_DST);
218 dev.queue.writeBuffer(this.bufAb, 0, new Float32Array(ab));
219 dev.queue.writeBuffer(this.bufAmm, 0, new Float32Array(amm));
220 dev.queue.writeBuffer(this.bufCtstw, 0, ctstw);
221 dev.queue.writeBuffer(this.bufTrig, 0, trig);
223 // --- shaders / pipelines ---
53ff7c5leg_analys: reduce with subgroupAdd, and retune wgAnalys for itdanfortunato 224 const subgroups = tuning('SHT_SUBGROUPS') !== false && dev.features.has('subgroups');
eb4a9e7Fold north/south latitude pairs onto one Legendre recurrencedanfortunato 225 // parity folding needs an equator-symmetric grid; Gauss nodes are, if nlat is even
226 const parity = tuning('SHT_PARITY') !== false && nlat % 2 === 0;
53ff7c5leg_analys: reduce with subgroupAdd, and retune wgAnalys for itdanfortunato 227 const wgAnalys =
228 (tuning('SHT_WG_ANALYS') as number | undefined) ??
229 defaultWgAnalys(nlat, dev.limits.maxComputeInvocationsPerWorkgroup, subgroups);
230 const legP = {
231 lmax,
232 mmax,
233 nlat,
234 wgSynth: WG_SYNTH,
235 wgAnalys,
236 subgroups,
b81424bleg_analys: reduce once per span of l-pairs, not once per pairdanfortunato 237 spanPairs: tuning('SHT_SPAN_PAIRS') as number | undefined,
eb4a9e7Fold north/south latitude pairs onto one Legendre recurrencedanfortunato 240 (this as { legLat: number }).legLat = parity ? nlat / 2 : nlat;
2dedc35turing-sphere: reaction-diffusion on the sphere, spectral solver on WebGPUJeremy Magland 241 const fourP = { mmax, nlat, nphi };
242 const [pLegS, pLegA, pFourS, pFourA] = await Promise.all([
243 makePipeline(dev, legSynthWGSL(legP), 'leg_synth'),
244 makePipeline(dev, legAnalysWGSL(legP), 'leg_analys'),
245 makePipeline(
246 dev,
247 this.fourierMode === 'fft' ? fftSynthWGSL(fourP) : dftSynthWGSL(fourP),
248 this.fourierMode === 'fft' ? 'fft_synth' : 'dft_synth',
249 ),
250 makePipeline(
251 dev,
252 this.fourierMode === 'fft' ? fftAnalysWGSL(fourP) : dftAnalysWGSL(fourP),
253 this.fourierMode === 'fft' ? 'fft_analys' : 'dft_analys',
254 ),
255 ]);
256 this.pipeLegSynth = pLegS;
257 this.pipeLegAnalys = pLegA;
258 this.pipeFourSynth = pFourS;
259 this.pipeFourAnalys = pFourA;
61e12f1Write the solver in MATLAB and compile it to WebGPUJeremy Magland 261 const entries = bgEntries;
2dedc35turing-sphere: reaction-diffusion on the sphere, spectral solver on WebGPUJeremy Magland 262 this.bgLegSynth = dev.createBindGroup({
263 layout: pLegS.getBindGroupLayout(0),
264 entries: entries([this.bufAb, this.bufAmm, this.bufCtstw, this.qlmIn, this.fmBuf]),
265 });
266 this.bgLegAnalys = dev.createBindGroup({
267 layout: pLegA.getBindGroupLayout(0),
268 entries: entries([this.bufAb, this.bufAmm, this.bufCtstw, this.fmBuf, this.qlmOut]),
269 });
270 this.bgFourSynth = dev.createBindGroup({
271 layout: pFourS.getBindGroupLayout(0),
272 entries: entries([this.fmBuf, this.spatBuf, this.bufTrig]),
273 });
274 this.bgFourAnalys = dev.createBindGroup({
275 layout: pFourA.getBindGroupLayout(0),
276 entries: entries([this.spatBuf, this.fmBuf, this.bufTrig]),
277 });
278 }
281 * Bind groups for one transform against caller-supplied spectral/spatial
282 * buffers, so a transform can read and write buffers it does not own (the
283 * .m-driven executor keeps a buffer per IR variable). Build these once at
284 * plan time, not per step. `fmBuf` stays internal scratch: passes and
285 * dispatches within a submission execute in order, so sequential transforms
286 * can share it.
287 */
288 createSynthBinding(qlmIn: GPUBuffer, spatOut: GPUBuffer): ShtBinding {
289 return {
290 bgLeg: this.device.createBindGroup({
291 layout: this.pipeLegSynth.getBindGroupLayout(0),
292 entries: bgEntries([this.bufAb, this.bufAmm, this.bufCtstw, qlmIn, this.fmBuf]),
293 }),
294 bgFour: this.device.createBindGroup({
295 layout: this.pipeFourSynth.getBindGroupLayout(0),
296 entries: bgEntries([this.fmBuf, spatOut, this.bufTrig]),
297 }),
298 };
299 }
301 createAnalysBinding(spatIn: GPUBuffer, qlmOut: GPUBuffer): ShtBinding {
302 return {
303 bgFour: this.device.createBindGroup({
304 layout: this.pipeFourAnalys.getBindGroupLayout(0),
305 entries: bgEntries([spatIn, this.fmBuf, this.bufTrig]),
306 }),
307 bgLeg: this.device.createBindGroup({
308 layout: this.pipeLegAnalys.getBindGroupLayout(0),
309 entries: bgEntries([this.bufAb, this.bufAmm, this.bufCtstw, this.fmBuf, qlmOut]),
310 }),
311 };
312 }
314 /** Record synthesis into an existing compute pass. */
315 encodeSynthInto(pass: GPUComputePassEncoder, b: ShtBinding): void {
2dedc35turing-sphere: reaction-diffusion on the sphere, spectral solver on WebGPUJeremy Magland 316 const { mmax, nlat, nphi } = this.cfg;
317 pass.setPipeline(this.pipeLegSynth);
61e12f1Write the solver in MATLAB and compile it to WebGPUJeremy Magland 318 pass.setBindGroup(0, b.bgLeg);
eb4a9e7Fold north/south latitude pairs onto one Legendre recurrencedanfortunato 319 pass.dispatchWorkgroups(Math.ceil(this.legLat / WG_SYNTH), mmax + 1);
2dedc35turing-sphere: reaction-diffusion on the sphere, spectral solver on WebGPUJeremy Magland 320 pass.setPipeline(this.pipeFourSynth);
61e12f1Write the solver in MATLAB and compile it to WebGPUJeremy Magland 321 pass.setBindGroup(0, b.bgFour);
2dedc35turing-sphere: reaction-diffusion on the sphere, spectral solver on WebGPUJeremy Magland 322 if (this.fourierMode === 'fft') {
323 pass.dispatchWorkgroups(nlat);
324 } else {
325 pass.dispatchWorkgroups(Math.ceil(nphi / 64), nlat);
326 }
327 }
61e12f1Write the solver in MATLAB and compile it to WebGPUJeremy Magland 329 /** Record analysis into an existing compute pass. */
330 encodeAnalysInto(pass: GPUComputePassEncoder, b: ShtBinding): void {
2dedc35turing-sphere: reaction-diffusion on the sphere, spectral solver on WebGPUJeremy Magland 331 const { mmax, nlat } = this.cfg;
332 pass.setPipeline(this.pipeFourAnalys);
61e12f1Write the solver in MATLAB and compile it to WebGPUJeremy Magland 333 pass.setBindGroup(0, b.bgFour);
2dedc35turing-sphere: reaction-diffusion on the sphere, spectral solver on WebGPUJeremy Magland 334 if (this.fourierMode === 'fft') {
335 pass.dispatchWorkgroups(nlat);
336 } else {
337 pass.dispatchWorkgroups(Math.ceil((mmax + 1) / 64), nlat);
338 }
339 pass.setPipeline(this.pipeLegAnalys);
61e12f1Write the solver in MATLAB and compile it to WebGPUJeremy Magland 340 pass.setBindGroup(0, b.bgLeg);
2dedc35turing-sphere: reaction-diffusion on the sphere, spectral solver on WebGPUJeremy Magland 341 pass.dispatchWorkgroups(mmax + 1);
344 /** Record the synthesis (spectral qlmIn -> spatial spatBuf) into an encoder. */
345 encodeSynth(encoder: GPUCommandEncoder): void {
346 const pass = encoder.beginComputePass({ label: 'sht-synth' });
347 this.encodeSynthInto(pass, { bgLeg: this.bgLegSynth, bgFour: this.bgFourSynth });
348 pass.end();
349 }
351 /** Record the analysis (spatial spatBuf -> spectral qlmOut) into an encoder. */
352 encodeAnalys(encoder: GPUCommandEncoder): void {
353 const pass = encoder.beginComputePass({ label: 'sht-analys' });
354 this.encodeAnalysInto(pass, { bgLeg: this.bgLegAnalys, bgFour: this.bgFourAnalys });
2dedc35turing-sphere: reaction-diffusion on the sphere, spectral solver on WebGPUJeremy Magland 355 pass.end();
356 }
358 /**
359 * Spectral -> spatial. qlm: interleaved [re, im], SHTNS LM ordering,
360 * length 2*nlm. Returns the spatial field, length nlat*nphi.
361 */
362 async synth(qlm: Float32Array): Promise<Float32Array> {
363 const { nlat, nphi } = this.cfg;
364 if (qlm.length !== 2 * this.nlm) throw new Error(`qlm must have length ${2 * this.nlm}`);
365 this.device.queue.writeBuffer(this.qlmIn, 0, qlm as Float32Array<ArrayBuffer>);
366 const enc = this.device.createCommandEncoder();
367 this.encodeSynth(enc);
368 enc.copyBufferToBuffer(this.spatBuf, 0, this.stageSpat, 0, 4 * nlat * nphi);
369 this.device.queue.submit([enc.finish()]);
370 await this.stageSpat.mapAsync(GPUMapMode.READ);
371 const out = new Float32Array(this.stageSpat.getMappedRange().slice(0));
372 this.stageSpat.unmap();
373 return out;
374 }
376 /** Spatial -> spectral. spat: length nlat*nphi. Returns interleaved qlm, length 2*nlm. */
377 async analys(spat: Float32Array): Promise<Float32Array> {
378 const { nlat, nphi } = this.cfg;
379 if (spat.length !== nlat * nphi) throw new Error(`spat must have length ${nlat * nphi}`);
380 this.device.queue.writeBuffer(this.spatBuf, 0, spat as Float32Array<ArrayBuffer>);
381 const enc = this.device.createCommandEncoder();
382 this.encodeAnalys(enc);
383 enc.copyBufferToBuffer(this.qlmOut, 0, this.stageQ, 0, 8 * this.nlm);
384 this.device.queue.submit([enc.finish()]);
385 await this.stageQ.mapAsync(GPUMapMode.READ);
386 const out = new Float32Array(this.stageQ.getMappedRange().slice(0));
387 this.stageQ.unmap();
388 return out;
389 }
391 destroy(): void {
392 for (const b of [
393 this.bufAb, this.bufAmm, this.bufCtstw, this.bufTrig, this.qlmIn, this.qlmOut,
394 this.fmBuf, this.spatBuf, this.stageSpat, this.stageQ,
395 ]) b?.destroy();
396 }
397}
35d91faDelete the TypeScript solver; the .m models are the only implementationJeremy Magland 399/** Best-effort human-readable adapter name, so it is clear which GPU (or
400 * software rasterizer) is actually running the transforms. */
401export async function describeAdapter(device: GPUDevice): Promise<string> {
402 const fmt = (info: GPUAdapterInfo | undefined): string => {
403 if (!info) return '';
404 const parts = [info.description, info.device, info.vendor].filter(
405 (s): s is string => !!s && s.length > 0,
406 );
407 const name = parts[0] ?? '';
408 return info.architecture && !name.includes(info.architecture)
409 ? `${name} (${info.architecture})`.trim()
410 : name;
411 };
412 const own = fmt((device as GPUDevice & { adapterInfo?: GPUAdapterInfo }).adapterInfo);
413 if (own) return own;
414 try {
415 const adapter = await navigator.gpu.requestAdapter();
416 return fmt(adapter?.info);
417 } catch {
418 return '';
419 }
420}
2dedc35turing-sphere: reaction-diffusion on the sphere, spectral solver on WebGPUJeremy Magland 422/** Request an adapter/device suitable for the transforms. */
423export async function requestShtDevice(): Promise<GPUDevice> {
424 if (!navigator.gpu) throw new Error('WebGPU is not available in this browser');
425 const adapter = await navigator.gpu.requestAdapter();
426 if (!adapter) throw new Error('No WebGPU adapter available');
427 // ask for a larger workgroup storage if the adapter offers it (bigger FFTs)
428 const wgStorage = Math.min(adapter.limits.maxComputeWorkgroupStorageSize, 32768);
53ff7c5leg_analys: reduce with subgroupAdd, and retune wgAnalys for itdanfortunato 429 // `subgroups` lets the analysis reduction use subgroupAdd instead of a
430 // shared-memory tree (2 barriers per l-pair instead of 1 + log2(wgAnalys)).
431 // Optional: ShtPlan falls back to the tree when it is not available.
432 const features: GPUFeatureName[] = adapter.features.has('subgroups') ? ['subgroups'] : [];
2dedc35turing-sphere: reaction-diffusion on the sphere, spectral solver on WebGPUJeremy Magland 433 return adapter.requestDevice({
53ff7c5leg_analys: reduce with subgroupAdd, and retune wgAnalys for itdanfortunato 434 requiredFeatures: features,
2dedc35turing-sphere: reaction-diffusion on the sphere, spectral solver on WebGPUJeremy Magland 435 requiredLimits: { maxComputeWorkgroupStorageSize: wgStorage },
436 });
437}