Sign in
concept-collection / fastandaccurate
Add mfs-gpu: the MFS on WebGPU, and a TypeScript form of the solver interface
The same method as mfs, with its assembly, its dense solve and its evaluation on a WebGPU device: a right-looking LU with partial pivoting, three dispatches per column, the right-hand side carried as an extra matrix column so that forward substitution disappears, about 3000 dispatches in one submit at n = 768. It runs on the page's own device in the browser and on Dawn from the command line, through the optional `webgpu` package; a run with no device skips it, the way a run with no matlab skips the MATLAB solvers. Since it cannot be a MATLAB file, the specification gains a second form of its solver interface (src/problems/laplace2d/problem.ts): the same curve, data and evaluation points as a plain object, under the same rules. Only the clock differs, a GPU solver having no tic/toc in its own runtime, so it is timed on the host around a run that ends by awaiting the device. Everything is f32, because WebGPU has no double precision, and for a conditioning-limited method that is the whole story: the error stops at 3.0e-7 on disk-easy where the same method in double reaches 9.5e-16, and the ceiling arrives at n = 96, where a dense solve is far too small to pay for a GPU. The curve is dominated everywhere by mfs. Emulating double in software, carrying each value as an unevaluated sum of two f32, was tried and does not survive this driver, which reassociates the error-free transformations that trick depends on. test/gpu-test.ts covers accuracy in node and scripts/check-gpu.mjs checks that the solver reaches a device from inside a real page and worker. Both skip where there is no device, which is what CI has.
Jeremy Magland <jmagland@flatironinstitute.org> committed commit acdea26b010e parent a69d634 Browse files
25 changed files+1415−91
README.mdmodified+29−12View file
@@ -13,15 +13,18 @@ problem-specified evaluation points. The central object is the
1313 resolution varies. No single ranking is presented; which curve wins can
1414 differ by accuracy regime, instance, and machine.
1515
16-Solvers are MATLAB function files. Most run via
16+Solvers are usually MATLAB function files. Most run via
1717 [numbl](https://numbl.org) (MATLAB syntax in the browser and in node),
1818 both on the site and from the command line; some run only in real
1919 MATLAB through the command line, and their results are marked as not
2020 reproducible in the browser. Two registry entries may share one file:
2121 the `-mat` solvers are their numbl twin's `solver.m` run in real MATLAB,
22-so that pair of curves measures the runtime rather than the method. Each
23-problem defines its own interface and instances in a written
24-specification; interfaces are per problem rather than shared.
22+so that pair of curves measures the runtime rather than the method. A
23+problem's interface also has a TypeScript form, for a solver that cannot
24+be a MATLAB file: `mfs-gpu` is the same method as `mfs` written in
25+TypeScript and WGSL and run on a WebGPU device. Each problem defines its
26+own interface and instances in a written specification; interfaces are
27+per problem rather than shared.
2528
2629 ## Problems
2730
@@ -60,6 +63,13 @@ makes its accelerated code path available without a Fortran compiler on
6063 the machine, since the mip fmm2d package ships a compiled MEX binary per
6164 platform.
6265
66+The solvers whose runtime is `webgpu` need a WebGPU device. In the
67+browser that is `navigator.gpu`; outside it, it is the optional
68+[`webgpu`](https://www.npmjs.com/package/webgpu) package (prebuilt Google
69+Dawn). That package is 68 MB, so the published command line does not ship
70+it and a run without it skips those solvers; `npm install webgpu` in a
71+checkout is enough to have them.
72+
6373 Useful flags: `--instance <id>`, `--solver <id>`, `--repeats N` (the
6474 minimum timed runs per point; each point is then repeated until it has
6575 used the `--time-budget`, 0.5 s by default), `--max-n N`, `--out dir`.
@@ -81,17 +91,24 @@ to this repository; see `src/solvers/`.
8191
8292 ```
8393 npm install
84-npm run dev # local dev server
85-npm test # solver convergence tests through numbl in node
86-npm run build # type-check, site build, CLI tarball (dist/)
87-npm run check-app # headless end-to-end check of the built site
94+npm run dev # local dev server
95+npm test # solver convergence tests through numbl in node
96+npm run test:matlab # the same for the MATLAB-runtime solvers (needs matlab)
97+npm run test:gpu # the same for the WebGPU solvers (needs a device)
98+npm run build # type-check, site build, CLI tarball (dist/)
99+npm run check-app # headless end-to-end check of the built site
100+npm run check-gpu # headless check that a WebGPU solver runs in a page
88101 ```
89102
103+Only `npm test` runs in CI, since a GitHub runner has neither MATLAB nor a
104+GPU; the other three skip cleanly where their runtime is missing and are
105+meant to be run locally before pushing solver changes.
106+
90107 Layout: `src/problems/` holds problem specs, instances, exact solutions,
91-and the problem-side MATLAB; `src/solvers/` the solver MATLAB files and
92-manifests; `src/harness/` the shared runner, sweep, and result schema
93-(used identically by the browser worker and the CLI); `src/app/` the
94-React site; `src/cli/` the command line.
108+the problem-side MATLAB and its TypeScript form; `src/solvers/` the
109+solver files and manifests; `src/harness/` the shared runner, sweep,
110+timing policy, and result schema (used identically by the browser worker
111+and the CLI); `src/app/` the React site; `src/cli/` the command line.
95112
96113 Deployed to GitHub Pages by `.github/workflows/deploy.yml` on push to
97114 main.
docs/problems/laplace-dirichlet-2d.mdmodified+37−0View file
@@ -223,6 +223,43 @@ The return value is a struct: `out.uEval` (m×1, required) and
223223 reconstruct the sources analytically or otherwise special-case the known
224224 solution; submissions are reviewed for this.
225225
226+### The TypeScript form
227+
228+A solver that cannot be a MATLAB file, one that runs on a GPU for
229+instance, receives the same information as a plain object instead, built by
230+`src/problems/laplace2d/problem.ts`:
231+
232+| field | meaning |
233+|---|---|
234+| `curve(t)` | `{x, y}`, the boundary point at parameter t |
235+| `curveD(t)`, `curveDD(t)` | its first and second derivatives |
236+| `g(t)` | Dirichlet data at boundary parameter t |
237+| `evalXY` | `Float64Array`, `nEval` points as interleaved x, y |
238+| `vizXY` | the visualization grid, interleaved, empty when not wanted |
239+
240+and returns `uEval`, and `uGrid` when the grid was asked for, as
241+`Float64Array`. Everything else is the same in both forms: the same
242+evaluation points in the same order, the same prohibition on
243+reconstructing the sources, and the same timing protocol. Only the clock
244+differs. A MATLAB solver is timed by tic/toc inside its own runtime; a GPU
245+solver has no such clock, so it is timed on the host around a run that
246+ends by awaiting the device, with the values read back before the clock
247+stops. That is the same synchronization point MATLAB's tic/toc gives.
248+Shader compilation and pipeline creation happen once per device rather
249+than once per run, which puts them where numbl's JIT compilation already
250+is: outside the timed runs, absorbed by the warmups.
251+
252+A note on what a GPU solver can and cannot compute here. WGSL has f32 and
253+f16 and no double precision, and no extension in the WebGPU standard adds
254+one. Emulating double in software (carrying a value as an unevaluated sum
255+of two f32) is not a reliable way out either: it rests on error-free
256+transformations such as `s = a + b; err = b - (s - a)`, which are exact
257+only if the compiler evaluates them as written, and WGSL permits an
258+implementation to use greater precision or to reassociate. So a WebGPU
259+solver on this problem works in single precision, and for a method whose
260+accuracy is limited by conditioning rather than by resolution that sets a
261+ceiling the CPU does not have. `mfs-gpu` measures where that ceiling is.
262+
226263 ## Visualization grid
227264
228265 When requested, `prob.vizXY` lists a 200×200 grid of points over the
package-lock.jsonmodified+25−1View file
@@ -11,16 +11,21 @@
1111 "fflate": "^0.8.3",
1212 "numbl": "^0.4.18",
1313 "react": "^19.2.7",
14- "react-dom": "^19.2.7"
14+ "react-dom": "^19.2.7",
15+ "webgpu": "^0.4.0"
1516 },
1617 "devDependencies": {
1718 "@types/node": "^24.0.0",
1819 "@types/react": "^19.2.0",
1920 "@types/react-dom": "^19.2.0",
2021 "@vitejs/plugin-react": "^5.0.0",
22+ "@webgpu/types": "^0.1.71",
2123 "tsx": "^4.19.0",
2224 "typescript": "~5.9.3",
2325 "vite": "^7.0.0"
26+ },
27+ "optionalDependencies": {
28+ "webgpu": "^0.4.0"
2429 }
2530 },
2631 "node_modules/@babel/code-frame": {
@@ -1383,6 +1388,13 @@
13831388 "vite": "^4.2.0 || ^5.0.0 || ^6.0.0 || ^7.0.0 || ^8.0.0"
13841389 }
13851390 },
1391+ "node_modules/@webgpu/types": {
1392+ "version": "0.1.71",
1393+ "resolved": "https://registry.npmjs.org/@webgpu/types/-/types-0.1.71.tgz",
1394+ "integrity": "sha512-mMy8/ODcKhab808co15eW+yN+HgXoQxRQHTiBV9Mrvl1r0ufnid7YOcI+gi4eUWSWl9ezD6TW2KXccrL8HCh2A==",
1395+ "devOptional": true,
1396+ "license": "BSD-3-Clause"
1397+ },
13861398 "node_modules/bail": {
13871399 "version": "2.0.2",
13881400 "resolved": "https://registry.npmjs.org/bail/-/bail-2.0.2.tgz",
@@ -3606,6 +3618,18 @@
36063618 }
36073619 }
36083620 },
3621+ "node_modules/webgpu": {
3622+ "version": "0.4.0",
3623+ "resolved": "https://registry.npmjs.org/webgpu/-/webgpu-0.4.0.tgz",
3624+ "integrity": "sha512-F5pimn3Aoi0zWjuRdiVs5TnrUwSzD2lESBohsIUsqyitWkGRQlXU2fhV6ycXlQTa1bvAf3sjqiUpBEpmSQ5ptA==",
3625+ "hasInstallScript": true,
3626+ "license": "MIT",
3627+ "optional": true,
3628+ "dependencies": {
3629+ "@webgpu/types": "^0.1.69",
3630+ "debug": "^4.4.0"
3631+ }
3632+ },
36093633 "node_modules/yallist": {
36103634 "version": "3.1.1",
36113635 "resolved": "https://registry.npmjs.org/yallist/-/yallist-3.1.1.tgz",
package.jsonmodified+7−1View file
@@ -10,7 +10,9 @@
1010 "preview": "vite preview",
1111 "test": "tsx test/solver-test.ts",
1212 "check-app": "node scripts/check-app.mjs",
13- "test:matlab": "tsx test/matlab-test.ts"
13+ "test:matlab": "tsx test/matlab-test.ts",
14+ "test:gpu": "tsx test/gpu-test.ts",
15+ "check-gpu": "node scripts/check-gpu.mjs"
1416 },
1517 "dependencies": {
1618 "fflate": "^0.8.3",
@@ -23,8 +25,12 @@
2325 "@types/react": "^19.2.0",
2426 "@types/react-dom": "^19.2.0",
2527 "@vitejs/plugin-react": "^5.0.0",
28+ "@webgpu/types": "^0.1.71",
2629 "tsx": "^4.19.0",
2730 "typescript": "~5.9.3",
2831 "vite": "^7.0.0"
32+ },
33+ "optionalDependencies": {
34+ "webgpu": "^0.4.0"
2935 }
3036 }
scripts/check-gpu.mjsadded+138−0View file
@@ -0,0 +1,138 @@
1+/**
2+ * End-to-end check of the WebGPU solver in a real browser: serves dist/,
3+ * opens the problem page in headless Chrome, picks mfs-gpu in the Solution
4+ * section, computes one point through the worker, and asserts the reported
5+ * error is what single precision gives. The node suite
6+ * (test/gpu-test.ts) covers accuracy across the instances; what this adds
7+ * is that the same code reaches a device from inside a page and a worker,
8+ * which is the only place the site's "Run in this browser" button lives.
9+ *
10+ * Usage: npm run build && node scripts/check-gpu.mjs
11+ * Exits 2, without failing, when no WebGPU device can be had here.
12+ */
13+import { createServer } from "node:http";
14+import { readFile } from "node:fs/promises";
15+import { join, extname } from "node:path";
16+import puppeteer from "puppeteer-core";
17+
18+const root = new URL("../dist", import.meta.url).pathname;
19+const types = {
20+ ".html": "text/html",
21+ ".js": "text/javascript",
22+ ".css": "text/css",
23+ ".wasm": "application/wasm",
24+ ".json": "application/json",
25+};
26+
27+const server = createServer(async (req, res) => {
28+ const path = req.url === "/" ? "/index.html" : (req.url ?? "/").split("?")[0];
29+ try {
30+ const data = await readFile(join(root, path));
31+ res.writeHead(200, {
32+ "content-type": types[extname(path)] ?? "application/octet-stream",
33+ });
34+ res.end(data);
35+ } catch {
36+ res.writeHead(404);
37+ res.end("not found");
38+ }
39+});
40+await new Promise((resolve) => server.listen(0, "127.0.0.1", resolve));
41+const port = server.address().port;
42+
43+// Hardware WebGPU first, then the software adapter, as turing-surface does.
44+const flagSets = [
45+ ["--headless=new", "--no-sandbox", "--enable-unsafe-webgpu", "--enable-features=Vulkan"],
46+ [
47+ "--headless=new",
48+ "--no-sandbox",
49+ "--enable-unsafe-webgpu",
50+ "--use-webgpu-adapter=swiftshader",
51+ "--enable-unsafe-swiftshader",
52+ ],
53+];
54+
55+/** f32 puts mfs-gpu's error on star-hard in this range; the point of the
56+ * check is that it computed at all, so the window is generous. */
57+const MIN_ERROR = 1e-9;
58+const MAX_ERROR = 1e-1;
59+
60+let outcome = null;
61+for (const flags of flagSets) {
62+ const browser = await puppeteer.launch({
63+ executablePath: process.env.CHROME_PATH ?? "/usr/bin/google-chrome",
64+ args: [...flags],
65+ protocolTimeout: 600_000,
66+ });
67+ try {
68+ const page = await browser.newPage();
69+ const pageErrors = [];
70+ page.on("pageerror", (err) => pageErrors.push(String(err)));
71+ await page.goto(`http://127.0.0.1:${port}/#/problem/laplace-dirichlet-2d`, {
72+ waitUntil: "networkidle2",
73+ });
74+ await page.waitForFunction(
75+ () => document.body.innerText.includes("Work-precision results"),
76+ { timeout: 60_000 }
77+ );
78+
79+ const gpu = await page.evaluate(() => "gpu" in navigator);
80+ if (!gpu) throw new Error("navigator.gpu absent");
81+
82+ // The Solution section's solver select, switched to mfs-gpu.
83+ const picked = await page.evaluate(() => {
84+ for (const sel of document.querySelectorAll("select")) {
85+ if ([...sel.options].some((o) => o.value === "mfs-gpu")) {
86+ sel.value = "mfs-gpu";
87+ sel.dispatchEvent(new Event("change", { bubbles: true }));
88+ return true;
89+ }
90+ }
91+ return false;
92+ });
93+ if (!picked) throw new Error("no solver select offering mfs-gpu");
94+
95+ const buttons = await page.$$("button");
96+ let compute = null;
97+ for (const b of buttons) {
98+ const t = await b.evaluate((el) => el.textContent);
99+ if (t && t.includes("Compute in this browser")) compute = b;
100+ }
101+ if (!compute) throw new Error("compute button not found");
102+ await compute.click();
103+
104+ const text = await page.waitForFunction(
105+ () => {
106+ const t = document.body.innerText;
107+ const m = /rel max error[^0-9eE.+-]*([0-9.]+e[+-][0-9]+)/i.exec(t);
108+ return m ? m[1] : false;
109+ },
110+ { timeout: 300_000 }
111+ );
112+ const relMax = Number(await text.jsonValue());
113+ if (pageErrors.length > 0) throw new Error(`page errors: ${pageErrors[0]}`);
114+ outcome = { flags, relMax };
115+ } catch (e) {
116+ console.error(`run with [${flags.join(" ")}] failed: ${e.message}`);
117+ } finally {
118+ await browser.close();
119+ }
120+ if (outcome) break;
121+ console.log("retrying with the software adapter…");
122+}
123+server.close();
124+
125+if (!outcome) {
126+ console.error("check-gpu: no WebGPU device in this browser; not run");
127+ process.exit(2);
128+}
129+if (!(outcome.relMax > MIN_ERROR && outcome.relMax < MAX_ERROR)) {
130+ console.error(
131+ `check-gpu: FAIL, mfs-gpu in the browser reported relMax ${outcome.relMax}, ` +
132+ `outside [${MIN_ERROR}, ${MAX_ERROR}]`
133+ );
134+ process.exit(1);
135+}
136+console.log(
137+ `check-gpu: mfs-gpu ran in the browser (relMax ${outcome.relMax.toExponential(2)})`
138+);
src/app/components/SolutionSection.tsxmodified+5−1View file
@@ -34,7 +34,11 @@ interface Computed {
3434 }
3535
3636 // Only solvers that run in the browser can compute a field here.
37-const BROWSER_SOLVERS = SOLVERS.filter((s) => s.runtime === "numbl");
37+// Everything a browser can run: numbl in the worker, and the WebGPU
38+// solvers on the worker's own device.
39+const BROWSER_SOLVERS = SOLVERS.filter(
40+ (s) => s.runtime === "numbl" || s.runtime === "webgpu"
41+);
3842
3943 export function SolutionSection({ inst }: { inst: Laplace2dInstance }) {
4044 const [solverId, setSolverId] = useState(BROWSER_SOLVERS[0].id);
src/app/matlabSources.tsmodified+30−10View file
@@ -1,18 +1,28 @@
1-// The MATLAB sources, inlined into the bundle by vite. Used by the worker
2-// and by the problem page's source listing; the node CLI reads the same
3-// files from disk.
1+// The solver and problem sources, inlined into the bundle by vite. Used by
2+// the worker (which needs the MATLAB text to run a numbl solver) and by the
3+// problem page's source listing, which also shows the TypeScript and WGSL
4+// of a WebGPU solver. The node CLI reads the same files from disk.
45
5-import { getSolver, solverSourceDir } from "../solvers";
6+import { getSolver, solverFiles, solverSourceDir } from "../solvers";
67
7-const files = import.meta.glob("../{problems,solvers}/**/*.m", {
8- query: "?raw",
9- import: "default",
10- eager: true,
11-}) as Record<string, string>;
8+const files = {
9+ ...(import.meta.glob("../{problems,solvers}/**/*.m", {
10+ query: "?raw",
11+ import: "default",
12+ eager: true,
13+ }) as Record<string, string>),
14+ // Only a WebGPU solver's own directory, not the registry modules that sit
15+ // one level up.
16+ ...(import.meta.glob("../solvers/*/*.ts", {
17+ query: "?raw",
18+ import: "default",
19+ eager: true,
20+ }) as Record<string, string>),
21+};
1222
1323 function get(path: string): string {
1424 const src = files[path];
15- if (!src) throw new Error(`missing MATLAB source: ${path}`);
25+ if (!src) throw new Error(`missing solver source: ${path}`);
1626 return src;
1727 }
1828
@@ -28,3 +38,13 @@ export function matlabBase() {
2838 export function solverSource(solverId: string): string {
2939 return get(`../solvers/${solverSourceDir(getSolver(solverId))}/solver.m`);
3040 }
41+
42+/** Every source file of a solver, named, for the page's listing. */
43+export function solverSourceFiles(
44+ solverId: string
45+): { name: string; code: string }[] {
46+ return solverFiles(getSolver(solverId)).map((f) => ({
47+ name: f.split("/").pop() as string,
48+ code: get(`../solvers/${f}`),
49+ }));
50+}
src/app/pages/AboutPage.tsxmodified+11−3View file
@@ -18,11 +18,14 @@ export function AboutPage() {
1818 points. A problem defines its own solver interface and a short list
1919 of official <strong>instances</strong> (parameter combinations) in a
2020 written specification, so every solver is compared on identical
21- inputs. Solvers are MATLAB function files. Most run via{" "}
21+ inputs. Solvers are usually MATLAB function files. Most run via{" "}
2222 <a href="https://numbl.org">numbl</a>, in the browser and from the
2323 command line alike; some run only in real MATLAB through the command
2424 line, and their results are marked as not reproducible in the
25- browser.
25+ browser. An interface also has a TypeScript form, for a solver that
26+ cannot be a MATLAB file: <code>mfs-gpu</code> is the same method as{" "}
27+ <code>mfs</code>, written in TypeScript and WGSL and run on a WebGPU
28+ device.
2629 </p>
2730
2831 <h2>Measurement</h2>
@@ -81,7 +84,12 @@ export function AboutPage() {
8184 <code>chunkie-dlp</code> needs one thing more, the{" "}
8285 <a href="https://mip.sh">mip</a> package manager on the MATLAB path,
8386 from which the harness installs chunkie and its FLAM and fmm2d
84- dependencies on first use.
87+ dependencies on first use. The solvers that run on WebGPU need a
88+ device: in the browser that is <code>navigator.gpu</code>, and
89+ outside it the optional{" "}
90+ <a href="https://www.npmjs.com/package/webgpu">webgpu</a> package
91+ (prebuilt Google Dawn), which is 68 MB and so is not shipped with
92+ the command line; a run without it skips them.
8593 </p>
8694 <p className="small muted">
8795 Note that npx caches by the exact URL string; the <code>?v=</code>{" "}
src/app/pages/ProblemPage.tsxmodified+25−18View file
@@ -15,7 +15,7 @@ import {
1515 import { solverColorVar } from "../colors";
1616 import { sweepInBrowser } from "../workerClient";
1717 import { DEFAULT_TIMING } from "../../harness/timing";
18-import { solverSource } from "../matlabSources";
18+import { solverSourceFiles } from "../matlabSources";
1919 import {
2020 WorkPrecisionChart,
2121 type ChartCurve,
@@ -265,12 +265,15 @@ export function ProblemPage({ problemId }: { problemId: string }) {
265265
266266 <h2>Solvers</h2>
267267 <p className="small muted" style={{ maxWidth: 640 }}>
268- Each solver is a MATLAB function file implementing the interface in
269- the <a href={SPEC_URL}>specification</a>. Most run through numbl, in
270- the browser and from the command line alike; the rest run only in
271- real MATLAB, through the command line. The <code>-mat</code> entries
272- are their numbl twin's file run that way, so each such pair measures
273- the runtime rather than the method.
268+ Most solvers are MATLAB function files implementing the interface in
269+ the <a href={SPEC_URL}>specification</a>, and most of those run
270+ through numbl, in the browser and from the command line alike; the
271+ rest run only in real MATLAB, through the command line. The{" "}
272+ <code>-mat</code> entries are their numbl twin's file run that way,
273+ so each such pair measures the runtime rather than the method.{" "}
274+ <code>mfs-gpu</code> is the same specification's second form:
275+ TypeScript and WGSL against the problem as a plain object, running on
276+ a WebGPU device.
274277 </p>
275278 {SOLVERS.map((s) => (
276279 <div
@@ -288,24 +291,28 @@ export function ProblemPage({ problemId }: { problemId: string }) {
288291 {s.id} v{s.version} · {s.backend} ·{" "}
289292 {s.runtime === "matlab"
290293 ? "runs in MATLAB via the command line"
291- : "runs via numbl in the browser and command line"}
294+ : s.runtime === "webgpu"
295+ ? "runs on WebGPU in the browser and command line"
296+ : "runs via numbl in the browser and command line"}
292297 {s.sourceDir ? ` · same solver.m as ${s.sourceDir}` : ""}
293298 </span>
294299 </div>
295300 <p className="small" style={{ color: "var(--text-2)" }}>
296301 {s.description}
297302 </p>
298- <details>
299- <summary className="small" style={{ cursor: "pointer" }}>
300- solver.m
301- </summary>
302- <pre style={{ maxHeight: 420, overflow: "auto", marginTop: 8 }}>
303- {solverSource(s.id)}
304- </pre>
305- </details>
303+ {solverSourceFiles(s.id).map((f) => (
304+ <details key={f.name}>
305+ <summary className="small" style={{ cursor: "pointer" }}>
306+ {f.name}
307+ </summary>
308+ <pre style={{ maxHeight: 420, overflow: "auto", marginTop: 8 }}>
309+ {f.code}
310+ </pre>
311+ </details>
312+ ))}
306313 <a
307314 className="small"
308- href={`${REPO_URL}/blob/main/src/solvers/${solverSourceDir(s)}/solver.m`}
315+ href={`${REPO_URL}/tree/main/src/solvers/${solverSourceDir(s)}`}
309316 >
310317 view on GitHub
311318 </a>
@@ -342,7 +349,7 @@ export function ProblemPage({ problemId }: { problemId: string }) {
342349 />
343350 {s.name}
344351 </label>{" "}
345- {s.runtime === "numbl" ? (
352+ {s.runtime === "numbl" || s.runtime === "webgpu" ? (
346353 <button
347354 onClick={() => runSolver(s.id)}
348355 disabled={running !== null}
src/app/results.tsmodified+14−4View file
@@ -48,10 +48,20 @@ export async function fetchCommittedResults(): Promise<ResultFile[]> {
4848 return results.filter((r): r is ResultFile => r !== null);
4949 }
5050
51-/** Short human label for the environment a result was measured in. */
51+/** Short human label for the environment a result was measured in. The
52+ * adapter is appended for a solver that ran on a GPU, since otherwise two
53+ * results from the same machine, one on its CPU and one on its GPU, would
54+ * carry the same label. */
5255 export function environmentLabel(r: ResultFile): string {
5356 const env = r.environment;
54- if (env.machineLabel) return `${env.machineLabel} (${env.kind})`;
55- if (env.kind === "browser") return "browser";
56- return `${env.cpu ?? "unknown cpu"} (${env.kind})`;
57+ const gpu = env.gpu ? ` on ${shortGpu(env.gpu)}` : "";
58+ if (env.machineLabel) return `${env.machineLabel} (${env.kind})${gpu}`;
59+ if (env.kind === "browser") return `browser${gpu}`;
60+ return `${env.cpu ?? "unknown cpu"} (${env.kind})${gpu}`;
61+}
62+
63+/** Adapter strings run long ("Intel open-source Mesa driver: Mesa 25.0.7"),
64+ * and a chart legend has no room for them. */
65+function shortGpu(s: string): string {
66+ return s.length > 28 ? `${s.slice(0, 27)}…` : s;
5767 }
src/app/worker.tsmodified+47−27View file
@@ -1,11 +1,14 @@
11 /// <reference lib="webworker" />
2-// The compute worker: runs numbl solves off the main thread. One request
3-// at a time; requests queue in the message queue while a sweep runs.
2+// The compute worker: runs solves off the main thread. One request at a
3+// time; requests queue in the message queue while a sweep runs. numbl
4+// solvers run synchronously here; a WebGPU solver runs on the worker's own
5+// device, which is why the handler is asynchronous.
46
57 import { getInstance } from "../problems/laplace2d/spec";
68 import { getSolver } from "../solvers";
79 import { runSweep } from "../harness/sweep";
8-import { runPoint } from "../harness/runner";
10+import { runPoint, type RunPoint } from "../harness/runner";
11+import { runPointGpu, runSweepGpu } from "../harness/webgpuRun";
912 import { toResultPoint, type ResultPoint } from "../harness/resultSchema";
1013 import { DEFAULT_TIMING, type TimingPolicy } from "../harness/timing";
1114 import { matlabBase, solverSource } from "./matlabSources";
@@ -39,26 +42,32 @@ export type WorkerResponse =
3942 }
4043 | { type: "error"; id: number; message: string };
4144
42-self.onmessage = (e: MessageEvent<WorkerRequest>) => {
45+self.onmessage = async (e: MessageEvent<WorkerRequest>) => {
4346 const msg = e.data;
4447 try {
4548 if (msg.type === "sweep") {
46- const points = runSweep({
47- instance: getInstance(msg.instanceId),
48- solver: getSolver(msg.solverId),
49- sources: { ...matlabBase(), solver: solverSource(msg.solverId) },
50- timing: msg.timing,
51- onPoint: (p, index, total) => {
52- const resp: WorkerResponse = {
53- type: "point",
54- id: msg.id,
55- point: toResultPoint(p),
56- index,
57- total,
58- };
59- postMessage(resp);
60- },
61- });
49+ const solver = getSolver(msg.solverId);
50+ const instance = getInstance(msg.instanceId);
51+ const onPoint = (p: RunPoint, index: number, total: number) => {
52+ const resp: WorkerResponse = {
53+ type: "point",
54+ id: msg.id,
55+ point: toResultPoint(p),
56+ index,
57+ total,
58+ };
59+ postMessage(resp);
60+ };
61+ const points =
62+ solver.runtime === "webgpu"
63+ ? await runSweepGpu({ instance, solver, timing: msg.timing, onPoint })
64+ : runSweep({
65+ instance,
66+ solver,
67+ sources: { ...matlabBase(), solver: solverSource(msg.solverId) },
68+ timing: msg.timing,
69+ onPoint,
70+ });
6271 const resp: WorkerResponse = {
6372 type: "sweepDone",
6473 id: msg.id,
@@ -66,13 +75,24 @@ self.onmessage = (e: MessageEvent<WorkerRequest>) => {
6675 };
6776 postMessage(resp);
6877 } else if (msg.type === "solution") {
69- const p = runPoint({
70- instance: getInstance(msg.instanceId),
71- n: msg.n,
72- timing: { ...DEFAULT_TIMING, minTimedRuns: 1, timeBudgetSeconds: 0 },
73- wantGrid: true,
74- sources: { ...matlabBase(), solver: solverSource(msg.solverId) },
75- });
78+ const solver = getSolver(msg.solverId);
79+ const timing = { ...DEFAULT_TIMING, minTimedRuns: 1, timeBudgetSeconds: 0 };
80+ const p =
81+ solver.runtime === "webgpu"
82+ ? await runPointGpu({
83+ instance: getInstance(msg.instanceId),
84+ solverId: msg.solverId,
85+ n: msg.n,
86+ timing,
87+ wantGrid: true,
88+ })
89+ : runPoint({
90+ instance: getInstance(msg.instanceId),
91+ n: msg.n,
92+ timing,
93+ wantGrid: true,
94+ sources: { ...matlabBase(), solver: solverSource(msg.solverId) },
95+ });
7696 if (!p.uGrid) throw new Error("solver returned no grid values");
7797 const resp: WorkerResponse = {
7898 type: "solutionDone",
src/cli/main.tsmodified+51−2View file
@@ -19,6 +19,9 @@ import { INSTANCES, getInstance } from "../problems/laplace2d/spec";
1919 import { setNumblFileIO } from "../harness/numblRun";
2020 import { NodeFileIOAdapter } from "./nodeFileIO";
2121 import { matlabAvailable, matlabSetup, runMatlabSweep } from "./matlabRun";
22+import { gpuUnavailableReason } from "../harness/webgpuDevice";
23+import { GPU_TIMER, runSweepGpu } from "../harness/webgpuRun";
24+import { getWebgpuSolver } from "../solvers/webgpuSolvers";
2225
2326 setNumblFileIO((vfs) => new NodeFileIOAdapter(vfs));
2427 import {
@@ -77,6 +80,23 @@ function environment(machineLabel: string | undefined, builtin: boolean): Result
7780 };
7881 }
7982
83+/** A WebGPU run records its adapter, and no numbl version: numbl is not
84+ * involved. */
85+function gpuEnvironment(
86+ machineLabel: string | undefined,
87+ gpu: string
88+): ResultEnvironment {
89+ return {
90+ kind: "node",
91+ runtime: `node ${process.version}`,
92+ gpu,
93+ os: `${os.platform()} ${os.release()}`,
94+ cpu: os.cpus()[0]?.model?.trim() ?? "unknown",
95+ machineLabel,
96+ browserReproducible: true,
97+ };
98+}
99+
80100 function matlabEnvironment(
81101 machineLabel: string | undefined,
82102 matlabVersion: string
@@ -190,11 +210,28 @@ async function runCommand(flags: Record<string, string>) {
190210 }
191211 wanted = wanted.filter((s) => s.runtime !== "matlab");
192212 }
213+ if (wanted.some((s) => s.runtime === "webgpu")) {
214+ const why = await gpuUnavailableReason();
215+ if (why !== null) {
216+ if (flags.solver) {
217+ throw new Error(`${flags.solver} runs on WebGPU: ${why}`);
218+ }
219+ for (const s of wanted.filter((x) => x.runtime === "webgpu")) {
220+ console.log(`skipping ${s.id}: runs on WebGPU. ${why}`);
221+ }
222+ wanted = wanted.filter((s) => s.runtime !== "webgpu");
223+ }
224+ }
193225 solverList = wanted.map((manifest) => ({
194226 manifest,
227+ // A WebGPU solver has no MATLAB source; the problem files are still
228+ // read so that the numbl and MATLAB entries share one code path.
195229 sources: {
196230 ...base,
197- solver: readSrc(`solvers/${solverSourceDir(manifest)}/solver.m`),
231+ solver:
232+ manifest.runtime === "webgpu"
233+ ? ""
234+ : readSrc(`solvers/${solverSourceDir(manifest)}/solver.m`),
198235 },
199236 source: "builtin",
200237 }));
@@ -216,7 +253,19 @@ async function runCommand(flags: Record<string, string>) {
216253 let resultPoints;
217254 let runEnv = env;
218255 let timer: string | undefined;
219- if (manifest.runtime === "matlab") {
256+ if (manifest.runtime === "webgpu") {
257+ const points = await runSweepGpu({
258+ instance: inst,
259+ solver: manifest,
260+ timing,
261+ maxN,
262+ onPoint: printPoint,
263+ });
264+ resultPoints = points.map(toResultPoint);
265+ const gpu = await getWebgpuSolver(manifest.id);
266+ runEnv = gpuEnvironment(flags.label, `${gpu.adapter} (${gpu.via})`);
267+ timer = GPU_TIMER;
268+ } else if (manifest.runtime === "matlab") {
220269 const setup = matlabSetup(manifest.id);
221270 const ns = sweepNFor(manifest, inst.id).filter(
222271 (n) => maxN === undefined || n <= maxN
src/harness/resultSchema.tsmodified+3−0View file
@@ -23,6 +23,9 @@ export interface ResultEnvironment {
2323 numblVersion?: string;
2424 os?: string;
2525 cpu?: string;
26+ /** The WebGPU adapter, for a solver that ran on one, and how WebGPU was
27+ * reached. Absent otherwise. */
28+ gpu?: string;
2629 /** Free-text label a human recognizes ("office workstation"). */
2730 machineLabel?: string;
2831 /** Whether a visitor can rerun this result in the browser. */
src/harness/webgpuDevice.tsadded+115−0View file
@@ -0,0 +1,115 @@
1+// Getting a WebGPU device, in the browser and outside it.
2+//
3+// In the browser navigator.gpu is there or it is not. In node it comes from
4+// the optional `webgpu` package (prebuilt Google Dawn), imported through a
5+// variable specifier so that neither the site bundle nor the command line
6+// bundle tries to resolve a native module at build time. The package is an
7+// optionalDependency and is 68 MB, so the command line does not ship it:
8+// a run without it skips the WebGPU solvers the same way a run without
9+// matlab on the PATH skips the MATLAB ones.
10+
11+export interface GpuEnvironment {
12+ device: GPUDevice;
13+ /** Adapter description for the result file's environment record. */
14+ adapter: string;
15+ /** How WebGPU was reached, for the same record. */
16+ via: string;
17+}
18+
19+let cached: Promise<GpuEnvironment> | null = null;
20+
21+/** Whether this is node rather than a page or a worker, which decides
22+ * whether a missing navigator.gpu means "install Dawn" or "this browser
23+ * does not have WebGPU". */
24+function isNode(): boolean {
25+ return typeof process !== "undefined" && !!process.versions?.node;
26+}
27+
28+async function installNodeWebGpu(): Promise<string> {
29+ const specifier = "webgpu";
30+ let mod: { create: (flags: string[]) => GPU; globals: Record<string, unknown> };
31+ try {
32+ mod = (await import(/* @vite-ignore */ specifier)) as typeof mod;
33+ } catch (e) {
34+ const detail = e instanceof Error ? e.message : String(e);
35+ if (/Cannot find (package|module) '?webgpu'?/.test(detail)) {
36+ throw new Error(
37+ "WebGPU outside the browser needs the optional `webgpu` package " +
38+ "(prebuilt Google Dawn): npm install webgpu"
39+ );
40+ }
41+ // Installed but unloadable is a different problem from missing, and
42+ // reporting it as missing sends people in circles.
43+ throw new Error(`the \`webgpu\` package is installed but did not load: ${detail}`);
44+ }
45+ Object.assign(globalThis, mod.globals);
46+ Object.defineProperty(globalThis, "navigator", {
47+ value: { gpu: mod.create([]) },
48+ configurable: true,
49+ writable: true,
50+ });
51+ return "node-webgpu (Google Dawn)";
52+}
53+
54+/** A device, requested once and shared. Throws with an actionable message
55+ * when there is no WebGPU here. */
56+export function requestGpu(): Promise<GpuEnvironment> {
57+ cached ??= (async () => {
58+ let via: string;
59+ if (typeof navigator !== "undefined" && navigator.gpu) {
60+ via = "browser";
61+ } else if (isNode()) {
62+ via = await installNodeWebGpu();
63+ } else {
64+ throw new Error(
65+ "this browser has no WebGPU: navigator.gpu is absent. Chrome and " +
66+ "Edge have it; Safari and Firefox need a recent version."
67+ );
68+ }
69+ const gpu = (navigator as Navigator).gpu;
70+ if (!gpu) throw new Error("no navigator.gpu after setup");
71+ const adapter = await gpu.requestAdapter();
72+ if (!adapter) {
73+ throw new Error(
74+ "WebGPU found no adapter. A headless machine often has none at all; " +
75+ "in Chrome chrome://gpu says why."
76+ );
77+ }
78+ const info = adapter.info as GPUAdapterInfo | undefined;
79+ const parts = [info?.vendor, info?.architecture, info?.device]
80+ .filter((s) => s)
81+ .join(" ");
82+ const device = await adapter.requestDevice();
83+ // A device lost mid-sweep would otherwise show up as a wrong answer.
84+ device.lost.then((reason) => {
85+ console.error(`WebGPU device lost: ${reason.reason} ${reason.message}`);
86+ });
87+ return {
88+ device,
89+ adapter: (info?.description || parts || "unknown adapter").trim(),
90+ via,
91+ };
92+ })();
93+ return cached;
94+}
95+
96+/** Whether a WebGPU device can be had here. Used to skip the WebGPU
97+ * solvers rather than fail a whole run. */
98+export async function gpuAvailable(): Promise<boolean> {
99+ try {
100+ await requestGpu();
101+ return true;
102+ } catch {
103+ return false;
104+ }
105+}
106+
107+/** Why WebGPU is unavailable, for a message to the user. */
108+export async function gpuUnavailableReason(): Promise<string | null> {
109+ try {
110+ await requestGpu();
111+ return null;
112+ } catch (e) {
113+ return e instanceof Error ? e.message : String(e);
114+ }
115+}
src/harness/webgpuRun.tsadded+99−0View file
@@ -0,0 +1,99 @@
1+// Runs a solver whose manifest declares runtime "webgpu": one that is
2+// TypeScript and WGSL rather than a MATLAB file, and executes on a WebGPU
3+// device. Used from the browser (where the device is the page's) and from
4+// the command line (where it is Dawn's, through the optional `webgpu`
5+// package).
6+//
7+// The protocol is the one in docs/problems/laplace-dirichlet-2d.md and the
8+// counting is identical to the numbl and MATLAB runners: one timed warmup
9+// reported as the cold time, a second untimed one, then timed runs until
10+// the policy is satisfied, of which the fastest is reported. What differs
11+// is the clock. A GPU solver has no tic/toc inside its own runtime, so the
12+// time is host wall clock around a run that ends by awaiting the device:
13+// the work is submitted and the result read back before the clock stops,
14+// which is the same synchronization point MATLAB's synchronous tic/toc
15+// gives. Shader compilation and pipeline creation happen once per device
16+// rather than per run, so, like numbl's JIT, they land outside the timed
17+// runs.
18+
19+import { buildProblem } from "../problems/laplace2d/problem";
20+import { evalErrors } from "../problems/laplace2d/exact";
21+import type { Laplace2dInstance } from "../problems/laplace2d/spec";
22+import { getWebgpuSolver } from "../solvers/webgpuSolvers";
23+import { sweepNFor, type SolverManifest } from "../solvers";
24+import { DEFAULT_TIMING, type TimingPolicy } from "./timing";
25+import type { RunPoint } from "./runner";
26+
27+export const GPU_TIMER = "host clock around submit and read-back";
28+
29+export interface GpuPointRequest {
30+ instance: Laplace2dInstance;
31+ solverId: string;
32+ n: number;
33+ timing?: TimingPolicy;
34+ wantGrid?: boolean;
35+}
36+
37+export async function runPointGpu(req: GpuPointRequest): Promise<RunPoint> {
38+ const timing = req.timing ?? DEFAULT_TIMING;
39+ const wantGrid = req.wantGrid ?? false;
40+ const solver = await getWebgpuSolver(req.solverId);
41+ const prob = buildProblem(req.instance, wantGrid);
42+ const { n } = req;
43+
44+ const t0 = performance.now();
45+ let out = await solver.run(prob, n, wantGrid);
46+ const coldSeconds = (performance.now() - t0) / 1000;
47+ out = await solver.run(prob, n, wantGrid);
48+
49+ const times: number[] = [];
50+ let total = 0;
51+ while (
52+ times.length < timing.maxTimedRuns &&
53+ (times.length < timing.minTimedRuns || total < timing.timeBudgetSeconds)
54+ ) {
55+ const t = performance.now();
56+ out = await solver.run(prob, n, wantGrid);
57+ const dt = (performance.now() - t) / 1000;
58+ times.push(dt);
59+ total += dt;
60+ }
61+
62+ const { relMax, relL2 } = evalErrors(req.instance, out.uEval);
63+ return {
64+ n,
65+ solveSeconds: Math.min(...times),
66+ solveSecondsAll: times,
67+ coldSeconds,
68+ relMax,
69+ relL2,
70+ uEval: out.uEval,
71+ uGrid: out.uGrid,
72+ };
73+}
74+
75+export interface GpuSweepOptions {
76+ instance: Laplace2dInstance;
77+ solver: SolverManifest;
78+ timing?: TimingPolicy;
79+ maxN?: number;
80+ onPoint?: (point: RunPoint, index: number, total: number) => void;
81+}
82+
83+export async function runSweepGpu(opts: GpuSweepOptions): Promise<RunPoint[]> {
84+ const ns = sweepNFor(opts.solver, opts.instance.id).filter(
85+ (n) => opts.maxN === undefined || n <= opts.maxN
86+ );
87+ const points: RunPoint[] = [];
88+ for (const [i, n] of ns.entries()) {
89+ const p = await runPointGpu({
90+ instance: opts.instance,
91+ solverId: opts.solver.id,
92+ n,
93+ timing: opts.timing,
94+ });
95+ points.push(p);
96+ opts.onPoint?.(p, i, ns.length);
97+ }
98+ return points;
99+}
src/problems/laplace2d/exact.tsmodified+19−0View file
@@ -37,6 +37,25 @@ export function boundaryRD(inst: Laplace2dInstance, t: number): number {
3737 return -inst.a * inst.k * Math.sin(inst.k * t);
3838 }
3939
40+/** d2r/dt2 of the boundary radius, written like boundaryRD through the
41+ * ratios f'/f and f''/f so that the rounded-square case stays away from
42+ * the underflow in f itself. */
43+export function boundaryRDD(inst: Laplace2dInstance, t: number): number {
44+ if (inst.shape === "rounded-square") {
45+ const p = inst.p as number;
46+ const c = Math.cos(t);
47+ const sn = Math.sin(t);
48+ const f = c ** p + sn ** p;
49+ const fp = p * (sn ** (p - 1) * c - c ** (p - 1) * sn);
50+ const fpp =
51+ p * ((p - 1) * (sn ** (p - 2) * c ** 2 + c ** (p - 2) * sn ** 2) - f);
52+ const u = fp / f;
53+ const v = fpp / f;
54+ return boundaryR(inst, t) * ((1 / p) * (1 / p + 1) * u ** 2 - v / p);
55+ }
56+ return -inst.a * inst.k * inst.k * Math.cos(inst.k * t);
57+}
58+
4059 /** The largest radius the boundary reaches, which sets the view extent
4160 * and the visualization grid. */
4261 export function maxRadius(inst: Laplace2dInstance): number {
src/problems/laplace2d/problem.tsadded+106−0View file
@@ -0,0 +1,106 @@
1+// The problem as a TypeScript object, for solvers that do not run through
2+// numbl. It carries exactly what build_problem.m hands a MATLAB solver —
3+// the curve with its derivatives, the Dirichlet data as a function of the
4+// boundary parameter, and the points where values are required — as plain
5+// functions and arrays. The specification
6+// (docs/problems/laplace-dirichlet-2d.md) states the interface once and
7+// this is its second form; a solver written against it, like the WebGPU
8+// MFS, sees the same information as a MATLAB one and no more. In
9+// particular the sources of the exact solution are used here only to
10+// manufacture g, and are not reachable from the returned object.
11+
12+import { type Laplace2dInstance } from "./spec";
13+import {
14+ boundaryR,
15+ boundaryRD,
16+ boundaryRDD,
17+ evalPoints,
18+ sources,
19+ vizGrid,
20+ VIZ_NGRID,
21+} from "./exact";
22+
23+export interface Vec2 {
24+ x: number;
25+ y: number;
26+}
27+
28+export interface Laplace2dProblem {
29+ /** Boundary point at parameter t. */
30+ curve(t: number): Vec2;
31+ /** First derivative of the curve with respect to t. */
32+ curveD(t: number): Vec2;
33+ /** Second derivative. */
34+ curveDD(t: number): Vec2;
35+ /** Dirichlet data at boundary parameter t. */
36+ g(t: number): number;
37+ /** The evaluation points: nEval rows of (x, y), interleaved. */
38+ evalXY: Float64Array;
39+ nEval: number;
40+ /** The visualization grid points, interleaved, empty when not wanted. */
41+ vizXY: Float64Array;
42+ nViz: number;
43+}
44+
45+/** Point and derivatives of x(t) = r(t) (cos t, sin t). */
46+function curveAt(inst: Laplace2dInstance, t: number, order: 0 | 1 | 2): Vec2 {
47+ const c = Math.cos(t);
48+ const s = Math.sin(t);
49+ const r = boundaryR(inst, t);
50+ if (order === 0) return { x: r * c, y: r * s };
51+ const r1 = boundaryRD(inst, t);
52+ if (order === 1) return { x: r1 * c - r * s, y: r1 * s + r * c };
53+ const r2 = boundaryRDD(inst, t);
54+ return {
55+ x: r2 * c - 2 * r1 * s - r * c,
56+ y: r2 * s + 2 * r1 * c - r * s,
57+ };
58+}
59+
60+export function buildProblem(
61+ inst: Laplace2dInstance,
62+ wantGrid = false
63+): Laplace2dProblem {
64+ // The three sources exist only to manufacture g, exactly as in
65+ // build_problem.m, and stay in this closure.
66+ const src = sources(inst);
67+ const pts = evalPoints(inst);
68+ const evalXY = new Float64Array(2 * pts.length);
69+ pts.forEach((p, i) => {
70+ evalXY[2 * i] = p.x;
71+ evalXY[2 * i + 1] = p.y;
72+ });
73+
74+ let vizXY = new Float64Array(0);
75+ if (wantGrid) {
76+ const { xs } = vizGrid(inst);
77+ vizXY = new Float64Array(2 * VIZ_NGRID * VIZ_NGRID);
78+ // Flat index p = ix * ngrid + iy, y varying fastest, matching
79+ // build_problem.m's meshgrid column order.
80+ let k = 0;
81+ for (let ix = 0; ix < VIZ_NGRID; ix++) {
82+ for (let iy = 0; iy < VIZ_NGRID; iy++) {
83+ vizXY[k++] = xs[ix];
84+ vizXY[k++] = xs[iy];
85+ }
86+ }
87+ }
88+
89+ return {
90+ curve: (t) => curveAt(inst, t, 0),
91+ curveD: (t) => curveAt(inst, t, 1),
92+ curveDD: (t) => curveAt(inst, t, 2),
93+ g: (t) => {
94+ const p = curveAt(inst, t, 0);
95+ let u = 0;
96+ for (const s of src) {
97+ u += s.c * 0.5 * Math.log((p.x - s.x) ** 2 + (p.y - s.y) ** 2);
98+ }
99+ return u;
100+ },
101+ evalXY,
102+ nEval: pts.length,
103+ vizXY,
104+ nViz: vizXY.length / 2,
105+ };
106+}
src/solvers/index.tsmodified+62−3View file
@@ -19,8 +19,11 @@ export interface SolverManifest {
1919 version: string;
2020 backend: "cpu" | "gpu";
2121 /** What executes the solver: "numbl" solvers run in the browser and in
22- * the CLI; "matlab" solvers run only in real MATLAB via the CLI. */
23- runtime: "numbl" | "matlab";
22+ * the CLI; "matlab" solvers run only in real MATLAB via the CLI;
23+ * "webgpu" solvers are TypeScript and WGSL rather than a MATLAB file and
24+ * run wherever a WebGPU device can be had, which is the browser and, with
25+ * the optional `webgpu` package, the CLI. */
26+ runtime: "numbl" | "matlab" | "webgpu";
2427 /** The resolution values a standard work-precision sweep runs. */
2528 sweepN: number[];
2629 /** Resolutions for instances that need a different range from sweepN,
@@ -34,11 +37,20 @@ export interface SolverManifest {
3437 sourceDir?: string;
3538 }
3639
37-/** Directory under src/solvers/ holding a solver's solver.m. */
40+/** Directory under src/solvers/ holding a solver's source. */
3841 export function solverSourceDir(s: SolverManifest): string {
3942 return s.sourceDir ?? s.id;
4043 }
4144
45+/** The solver's source files, relative to src/solvers/. A MATLAB solver is
46+ * one file; a WebGPU one is its TypeScript driver and the module that
47+ * generates its WGSL. */
48+export function solverFiles(s: SolverManifest): string[] {
49+ const dir = solverSourceDir(s);
50+ if (s.runtime === "webgpu") return [`${dir}/solver.ts`, `${dir}/wgsl.ts`];
51+ return [`${dir}/solver.m`];
52+}
53+
4254 /** The resolutions a sweep of this solver runs on this instance. */
4355 export function sweepNFor(s: SolverManifest, instanceId: string): number[] {
4456 return s.sweepNByInstance?.[instanceId] ?? s.sweepN;
@@ -85,6 +97,53 @@ export const SOLVERS: SolverManifest[] = [
8597 runtime: "numbl",
8698 sweepN: [8, 12, 16, 24, 32, 48, 64, 96, 128, 192, 256, 384, 512, 768],
8799 },
100+ {
101+ id: "mfs-gpu",
102+ name: "Method of fundamental solutions (WebGPU)",
103+ description:
104+ "The mfs method unchanged, with its assembly, its dense solve and " +
105+ "its evaluation all on the GPU through WebGPU. It is the one solver " +
106+ "here that is not a MATLAB file: TypeScript and WGSL, run on the " +
107+ "page's own device in the browser and on Dawn from the command line. " +
108+ "The solve is a right-looking LU with partial pivoting, three " +
109+ "dispatches per column, with the right-hand side carried as an extra " +
110+ "matrix column so that forward substitution disappears; at n = 768 " +
111+ "that is some 3000 dispatches in a single submit. Everything is in " +
112+ "f32, because WebGPU has no double precision and no extension in the " +
113+ "standard adds one. For a method whose accuracy is set by " +
114+ "conditioning rather than by resolution that is decisive: on " +
115+ "disk-easy the error stops at 3.0e-7 where the same method in double " +
116+ "reaches 9.5e-16, and likewise 4.0e-7 against 7.7e-16 on star-medium " +
117+ "and 2.7e-6 against 2.3e-13 on square-corners. On the three " +
118+ "instances where mfs is already conditioning-limited in double the " +
119+ "gap narrows to two and a half or three orders (1.4e-3 against " +
120+ "4.0e-6 on star-hard), since there both are losing to the same " +
121+ "ill-conditioning and only the rate differs. Note where the ceiling " +
122+ "arrives: n = 96 charges, which is far too small a dense solve to " +
123+ "pay for a GPU. At that resolution this solver is 4.5 times slower " +
124+ "than mfs in the same browser and ten times slower than mfs-mat in " +
125+ "real MATLAB, and the dispatch overhead is plain at the bottom of " +
126+ "the sweep, where 8 charges still cost a millisecond. It does become " +
127+ "the faster of the two browser solvers past about n = 384, reaching " +
128+ "6.6 times at n = 768, but f32 has taken its own error to 1e-5 by " +
129+ "then. So the curve is dominated everywhere: at every accuracy this " +
130+ "solver reaches, mfs on the CPU reaches it sooner. We emphasize that " +
131+ "this is a statement about the method and the size of its systems " +
132+ "rather than about the hardware. An O(n^3) solve does eventually pay " +
133+ "for a GPU, and the same kernels would; the MFS simply cannot use an " +
134+ "n that large, because single precision has ended its convergence " +
135+ "long before. A method that stays accurate as n grows, which on this " +
136+ "problem means the integral-equation solvers, is where a WebGPU " +
137+ "backend would have something to win.",
138+ version: "1.0.0",
139+ backend: "gpu",
140+ runtime: "webgpu",
141+ // The same list as mfs, so the two curves land on the same resolutions
142+ // and the pair isolates the backend. Most of the upper half is past
143+ // where f32 stops improving, and is there to show both the ceiling and
144+ // what the O(n^3) solve costs on either side.
145+ sweepN: [8, 12, 16, 24, 32, 48, 64, 96, 128, 192, 256, 384, 512, 768],
146+ },
88147 {
89148 id: "nystrom-dlp",
90149 name: "Nystrom double-layer BIE",
src/solvers/mfs-gpu/solver.tsadded+263−0View file
@@ -0,0 +1,263 @@
1+// Method of fundamental solutions on WebGPU.
2+//
3+// The same method as src/solvers/mfs/solver.m: n logarithmic point charges
4+// on a curve a fixed distance 0.3 outside the boundary, their strengths
5+// found by collocating the Dirichlet data at n boundary points, and the
6+// potential evaluated at the requested targets. What differs is where the
7+// work happens. The assembly, the dense solve, and the evaluation all run
8+// as WebGPU compute passes; the host does only the O(n) geometry, in f64,
9+// the way shtns-webgpu and the rest of the WebGPU work in this collection
10+// keep precomputation on the CPU.
11+//
12+// The solve is a right-looking LU with partial pivoting, three dispatches
13+// per column: pivot search and row swap in one workgroup, the multipliers,
14+// then the rank-one update of the trailing submatrix. There is no blocking
15+// and no GEMM, so it is memory-bound rather than compute-bound; it is the
16+// straightforward implementation, and a blocked panel factorization would
17+// be the next thing to try. The right-hand side rides along as an extra
18+// matrix column, which is what makes forward substitution disappear; back
19+// substitution is one dispatch per column after that. A sweep at n = 768
20+// therefore encodes about 3000 dispatches, all into one command buffer and
21+// one submit.
22+//
23+// Everything is f32: WebGPU has no double precision, which is what caps
24+// this solver's accuracy well short of the same method on the CPU. See
25+// ./wgsl.ts.
26+
27+import { requestGpu } from "../../harness/webgpuDevice";
28+import type { Laplace2dProblem } from "../../problems/laplace2d/problem";
29+import { EVAL_GROUP, LANES, mfsShader, TILE } from "./wgsl";
30+
31+/** Distance of the charge curve outside the boundary, as in mfs/solver.m. */
32+const DELTA = 0.3;
33+
34+const KERNELS = [
35+ "assemble",
36+ "pivot",
37+ "multipliers",
38+ "update",
39+ "backsub",
40+ "evaluate",
41+] as const;
42+type Kernel = (typeof KERNELS)[number];
43+
44+export interface GpuMfsResult {
45+ uEval: Float64Array;
46+ uGrid: Float64Array | null;
47+}
48+
49+function ceilDiv(a: number, b: number): number {
50+ return Math.ceil(a / b);
51+}
52+
53+export class MfsGpu {
54+ private constructor(
55+ private readonly device: GPUDevice,
56+ private readonly layout: GPUBindGroupLayout,
57+ private readonly pipelines: Record<Kernel, GPUComputePipeline>,
58+ private readonly stepStride: number,
59+ /** The adapter, for the result file's environment record. */
60+ readonly adapter: string,
61+ readonly via: string
62+ ) {}
63+
64+ /** Compile the shader and build the pipelines. Done once per device; a
65+ * solver call reuses them, which is where numbl's JIT sits too. */
66+ static async create(): Promise<MfsGpu> {
67+ const { device, adapter, via } = await requestGpu();
68+ const module = device.createShaderModule({
69+ code: mfsShader(),
70+ label: "mfs-gpu",
71+ });
72+ const info = await module.getCompilationInfo();
73+ const errors = info.messages.filter((m) => m.type === "error");
74+ if (errors.length > 0) {
75+ throw new Error(
76+ "mfs-gpu failed to compile:\n" +
77+ errors.map((m) => ` ${m.lineNum}:${m.linePos} ${m.message}`).join("\n")
78+ );
79+ }
80+ const layout = device.createBindGroupLayout({
81+ entries: [
82+ { binding: 0, visibility: GPUShaderStage.COMPUTE, buffer: { type: "uniform" } },
83+ {
84+ binding: 1,
85+ visibility: GPUShaderStage.COMPUTE,
86+ buffer: { type: "uniform", hasDynamicOffset: true, minBindingSize: 16 },
87+ },
88+ { binding: 2, visibility: GPUShaderStage.COMPUTE, buffer: { type: "read-only-storage" } },
89+ { binding: 3, visibility: GPUShaderStage.COMPUTE, buffer: { type: "storage" } },
90+ { binding: 4, visibility: GPUShaderStage.COMPUTE, buffer: { type: "storage" } },
91+ { binding: 5, visibility: GPUShaderStage.COMPUTE, buffer: { type: "storage" } },
92+ { binding: 6, visibility: GPUShaderStage.COMPUTE, buffer: { type: "storage" } },
93+ ],
94+ });
95+ const pipelineLayout = device.createPipelineLayout({ bindGroupLayouts: [layout] });
96+ const built = await Promise.all(
97+ KERNELS.map((k) =>
98+ device.createComputePipelineAsync({
99+ layout: pipelineLayout,
100+ compute: { module, entryPoint: k },
101+ label: `mfs-gpu ${k}`,
102+ })
103+ )
104+ );
105+ const pipelines = Object.fromEntries(
106+ KERNELS.map((k, i) => [k, built[i]])
107+ ) as Record<Kernel, GPUComputePipeline>;
108+ return new MfsGpu(
109+ device,
110+ layout,
111+ pipelines,
112+ Math.max(16, device.limits.minUniformBufferOffsetAlignment),
113+ adapter,
114+ via
115+ );
116+ }
117+
118+ private checked = false;
119+
120+ /**
121+ * One full solve at resolution n: the timed unit of the protocol, so it
122+ * includes the host geometry, the buffer allocation, the assembly, the
123+ * factorization, the evaluation, and the read-back.
124+ */
125+ async run(
126+ prob: Laplace2dProblem,
127+ n: number,
128+ wantGrid = false
129+ ): Promise<GpuMfsResult> {
130+ const { device } = this;
131+ const ld = n + 1;
132+ const nGrid = wantGrid ? prob.nViz : 0;
133+ const m = prob.nEval + nGrid;
134+ const bytes = 4;
135+
136+ // --- host geometry, in f64, rounded once on the way to the device --
137+ const data = new Float32Array(5 * n + 2 * m);
138+ const put = (index: number, v: number) => {
139+ data[index] = v;
140+ };
141+ for (let j = 0; j < n; j++) {
142+ const t = (2 * Math.PI * j) / n;
143+ const p = prob.curve(t);
144+ const d = prob.curveD(t);
145+ const sp = Math.hypot(d.x, d.y);
146+ // The outward unit normal of the counterclockwise curve.
147+ put(2 * j, p.x);
148+ put(2 * j + 1, p.y);
149+ put(2 * n + 2 * j, p.x + (DELTA * d.y) / sp);
150+ put(2 * n + 2 * j + 1, p.y - (DELTA * d.x) / sp);
151+ put(4 * n + j, prob.g(t));
152+ }
153+ for (let i = 0; i < prob.nEval; i++) {
154+ put(5 * n + 2 * i, prob.evalXY[2 * i]);
155+ put(5 * n + 2 * i + 1, prob.evalXY[2 * i + 1]);
156+ }
157+ for (let i = 0; i < nGrid; i++) {
158+ const o = 5 * n + 2 * (prob.nEval + i);
159+ put(o, prob.vizXY[2 * i]);
160+ put(o + 1, prob.vizXY[2 * i + 1]);
161+ }
162+
163+ // --- buffers -----------------------------------------------------
164+ const S = GPUBufferUsage.STORAGE;
165+ const buffers = {
166+ dims: device.createBuffer({ size: 16, usage: GPUBufferUsage.UNIFORM | GPUBufferUsage.COPY_DST }),
167+ step: device.createBuffer({
168+ size: Math.max(this.stepStride * Math.max(n, 1), this.stepStride),
169+ usage: GPUBufferUsage.UNIFORM | GPUBufferUsage.COPY_DST,
170+ }),
171+ data: device.createBuffer({ size: data.byteLength, usage: S | GPUBufferUsage.COPY_DST }),
172+ mat: device.createBuffer({ size: n * ld * bytes, usage: S }),
173+ lcol: device.createBuffer({ size: n * bytes, usage: S }),
174+ sol: device.createBuffer({ size: n * bytes, usage: S }),
175+ out: device.createBuffer({ size: m * bytes, usage: S | GPUBufferUsage.COPY_SRC }),
176+ read: device.createBuffer({
177+ size: m * bytes,
178+ usage: GPUBufferUsage.COPY_DST | GPUBufferUsage.MAP_READ,
179+ }),
180+ };
181+ try {
182+ device.queue.writeBuffer(buffers.dims, 0, new Uint32Array([n, ld, m, 0]));
183+ // One slot per elimination step, so k reaches the shader through a
184+ // dynamic uniform offset and the whole sweep needs one bind group.
185+ const steps = new Uint32Array((this.stepStride / 4) * Math.max(n, 1));
186+ for (let k = 0; k < n; k++) steps[k * (this.stepStride / 4)] = k;
187+ device.queue.writeBuffer(buffers.step, 0, steps);
188+ device.queue.writeBuffer(buffers.data, 0, data);
189+
190+ const bind = device.createBindGroup({
191+ layout: this.layout,
192+ entries: [
193+ { binding: 0, resource: { buffer: buffers.dims } },
194+ { binding: 1, resource: { buffer: buffers.step, size: 16 } },
195+ { binding: 2, resource: { buffer: buffers.data } },
196+ { binding: 3, resource: { buffer: buffers.mat } },
197+ { binding: 4, resource: { buffer: buffers.lcol } },
198+ { binding: 5, resource: { buffer: buffers.sol } },
199+ { binding: 6, resource: { buffer: buffers.out } },
200+ ],
201+ });
202+
203+ // The first run of a session is checked for validation errors, which
204+ // is a host round trip and so does not belong in a timed run; the
205+ // protocol's untimed warmups absorb it.
206+ const check = !this.checked;
207+ if (check) device.pushErrorScope("validation");
208+
209+ const enc = device.createCommandEncoder();
210+ const pass = enc.beginComputePass();
211+ const P = this.pipelines;
212+ // WebGPU orders dispatches within a pass and makes each one's writes
213+ // visible to the next, so the dependences here need no barriers.
214+ pass.setPipeline(P.assemble);
215+ pass.setBindGroup(0, bind, [0]);
216+ pass.dispatchWorkgroups(ceilDiv(ld, TILE), ceilDiv(n, TILE));
217+ for (let k = 0; k < n; k++) {
218+ const off = k * this.stepStride;
219+ pass.setPipeline(P.pivot);
220+ pass.setBindGroup(0, bind, [off]);
221+ pass.dispatchWorkgroups(1);
222+ const rows = n - 1 - k;
223+ if (rows > 0) {
224+ pass.setPipeline(P.multipliers);
225+ pass.setBindGroup(0, bind, [off]);
226+ pass.dispatchWorkgroups(ceilDiv(rows, LANES));
227+ pass.setPipeline(P.update);
228+ pass.setBindGroup(0, bind, [off]);
229+ pass.dispatchWorkgroups(ceilDiv(ld - 1 - k, TILE), ceilDiv(rows, TILE));
230+ }
231+ }
232+ for (let k = n - 1; k >= 0; k--) {
233+ pass.setPipeline(P.backsub);
234+ pass.setBindGroup(0, bind, [k * this.stepStride]);
235+ pass.dispatchWorkgroups(1);
236+ }
237+ pass.setPipeline(P.evaluate);
238+ pass.setBindGroup(0, bind, [0]);
239+ pass.dispatchWorkgroups(ceilDiv(m, EVAL_GROUP));
240+ pass.end();
241+ enc.copyBufferToBuffer(buffers.out, 0, buffers.read, 0, m * bytes);
242+ device.queue.submit([enc.finish()]);
243+
244+ await buffers.read.mapAsync(GPUMapMode.READ);
245+ const raw = new Float32Array(buffers.read.getMappedRange().slice(0));
246+ buffers.read.unmap();
247+
248+ if (check) {
249+ this.checked = true;
250+ const err = await device.popErrorScope();
251+ if (err) throw new Error(`mfs-gpu: WebGPU validation: ${err.message}`);
252+ }
253+
254+ const all = Float64Array.from(raw.subarray(0, m));
255+ return {
256+ uEval: all.slice(0, prob.nEval),
257+ uGrid: wantGrid ? all.slice(prob.nEval) : null,
258+ };
259+ } finally {
260+ for (const b of Object.values(buffers)) b.destroy();
261+ }
262+ }
263+}
src/solvers/mfs-gpu/wgsl.tsadded+177−0View file
@@ -0,0 +1,177 @@
1+// The WGSL for the WebGPU method of fundamental solutions.
2+//
3+// Everything here is f32, because that is what WebGPU has: WGSL's floating
4+// point types are f32 and f16, there is no f64, and no extension in the
5+// standard adds one. That single fact is the most interesting thing about
6+// this solver, because the MFS is conditioning-limited: its accuracy stops
7+// improving where rounding in the collocation solve overtakes the
8+// approximation error, and in f32 that happens six to nine orders short of
9+// where it happens in double. See the solver's manifest entry for the
10+// measured numbers.
11+//
12+// Note that software emulation of double precision (carrying each value as
13+// an unevaluated sum of two f32, the "double-single" trick) is not a way
14+// out. It depends on error-free transformations such as
15+// s = a + b; err = b - (s - a), which are only error-free if the compiler
16+// evaluates them exactly as written; WGSL permits an implementation to
17+// compute with greater precision or to reassociate, and the Mesa driver
18+// this was tried on does exactly that, returning f32-accurate results for
19+// values loaded from a buffer. So f32 is what a WebGPU solver gets.
20+
21+/** Workgroup sizes: the two-dimensional kernels (assembly and the trailing
22+ * update), the one-dimensional ones, and the evaluation. */
23+export const TILE = 16;
24+export const LANES = 256;
25+export const EVAL_GROUP = 64;
26+
27+/**
28+ * The whole solver as one shader module.
29+ *
30+ * The matrix is the augmented system [A | b], n rows of ld = n + 1 entries,
31+ * row-major, so that the elimination applies to the right-hand side as it
32+ * goes and no separate forward substitution is needed. The multipliers are
33+ * never stored: nothing reads them again.
34+ *
35+ * One read-only buffer holds every input, because the default WebGPU limit
36+ * is eight storage buffers per stage. Its layout, in f32 elements, is
37+ * derived from n rather than passed: collocation points (x, y interleaved)
38+ * at 0, charges at 2n, the Dirichlet data at 4n, the targets (interleaved)
39+ * at 5n.
40+ */
41+export function mfsShader(): string {
42+ return `// generated by src/solvers/mfs-gpu/wgsl.ts
43+
44+struct Dims {
45+ n: u32, // charges = collocation points
46+ ld: u32, // row stride of the augmented matrix, n + 1
47+ m: u32, // evaluation targets
48+ pad: u32,
49+};
50+
51+struct Step { k: u32, pad0: u32, pad1: u32, pad2: u32 };
52+
53+@group(0) @binding(0) var<uniform> dims: Dims;
54+@group(0) @binding(1) var<uniform> step: Step;
55+@group(0) @binding(2) var<storage, read> data: array<f32>;
56+@group(0) @binding(3) var<storage, read_write> mat: array<f32>;
57+@group(0) @binding(4) var<storage, read_write> lcol: array<f32>;
58+@group(0) @binding(5) var<storage, read_write> sol: array<f32>;
59+@group(0) @binding(6) var<storage, read_write> out: array<f32>;
60+
61+fn pt_x(i: u32) -> f32 { return data[2u * i]; }
62+fn pt_y(i: u32) -> f32 { return data[2u * i + 1u]; }
63+fn ch_x(j: u32) -> f32 { return data[2u * dims.n + 2u * j]; }
64+fn ch_y(j: u32) -> f32 { return data[2u * dims.n + 2u * j + 1u]; }
65+fn bdata(i: u32) -> f32 { return data[4u * dims.n + i]; }
66+fn tg_x(t: u32) -> f32 { return data[5u * dims.n + 2u * t]; }
67+fn tg_y(t: u32) -> f32 { return data[5u * dims.n + 2u * t + 1u]; }
68+
69+/** The MFS kernel log|p - q|, written as half the log of the squared
70+ * distance so that no square root is taken. */
71+fn logdist(ax: f32, ay: f32, bx: f32, by: f32) -> f32 {
72+ let dx = ax - bx;
73+ let dy = ay - by;
74+ return 0.5 * log(dx * dx + dy * dy);
75+}
76+
77+@compute @workgroup_size(${TILE}, ${TILE})
78+fn assemble(@builtin(global_invocation_id) gid: vec3<u32>) {
79+ let i = gid.y;
80+ let j = gid.x;
81+ if (i >= dims.n || j >= dims.ld) { return; }
82+ if (j == dims.n) {
83+ mat[i * dims.ld + j] = bdata(i);
84+ return;
85+ }
86+ mat[i * dims.ld + j] = logdist(pt_x(i), pt_y(i), ch_x(j), ch_y(j));
87+}
88+
89+var<workgroup> best_val: array<f32, ${LANES}>;
90+var<workgroup> best_idx: array<u32, ${LANES}>;
91+
92+/** Partial pivoting for step k: find the largest |A(i,k)| over i >= k and
93+ * swap that row with row k. One workgroup, so the search reduces in shared
94+ * memory and the swap needs no second dispatch. */
95+@compute @workgroup_size(${LANES})
96+fn pivot(@builtin(local_invocation_id) lid: vec3<u32>) {
97+ let k = step.k;
98+ let n = dims.n;
99+ let ld = dims.ld;
100+ var bv = -1.0;
101+ var bi = k;
102+ for (var i = k + lid.x; i < n; i = i + ${LANES}u) {
103+ let v = abs(mat[i * ld + k]);
104+ if (v > bv) { bv = v; bi = i; }
105+ }
106+ best_val[lid.x] = bv;
107+ best_idx[lid.x] = bi;
108+ workgroupBarrier();
109+ for (var s = ${LANES / 2}u; s > 0u; s = s >> 1u) {
110+ if (lid.x < s && best_val[lid.x + s] > best_val[lid.x]) {
111+ best_val[lid.x] = best_val[lid.x + s];
112+ best_idx[lid.x] = best_idx[lid.x + s];
113+ }
114+ workgroupBarrier();
115+ }
116+ let p = best_idx[0];
117+ if (p != k) {
118+ for (var j = k + lid.x; j < ld; j = j + ${LANES}u) {
119+ let a = mat[k * ld + j];
120+ mat[k * ld + j] = mat[p * ld + j];
121+ mat[p * ld + j] = a;
122+ }
123+ }
124+}
125+
126+/** The multipliers of step k, computed once rather than once per element of
127+ * the trailing update. */
128+@compute @workgroup_size(${LANES})
129+fn multipliers(@builtin(global_invocation_id) gid: vec3<u32>) {
130+ let k = step.k;
131+ let i = k + 1u + gid.x;
132+ if (i >= dims.n) { return; }
133+ lcol[i] = mat[i * dims.ld + k] / mat[k * dims.ld + k];
134+}
135+
136+/** The rank-one update of the trailing submatrix, columns k+1 .. ld-1, so
137+ * the augmented right-hand column is eliminated along with the rest. */
138+@compute @workgroup_size(${TILE}, ${TILE})
139+fn update(@builtin(global_invocation_id) gid: vec3<u32>) {
140+ let k = step.k;
141+ let i = k + 1u + gid.y;
142+ let j = k + 1u + gid.x;
143+ if (i >= dims.n || j >= dims.ld) { return; }
144+ let ld = dims.ld;
145+ mat[i * ld + j] = mat[i * ld + j] - lcol[i] * mat[k * ld + j];
146+}
147+
148+/** Back substitution, one step per dispatch: x(k) from row k of U and the
149+ * running right-hand column, then that column updated above row k. Every
150+ * lane recomputes x(k) so that no barrier is needed. */
151+@compute @workgroup_size(${LANES})
152+fn backsub(@builtin(local_invocation_id) lid: vec3<u32>) {
153+ let k = step.k;
154+ let ld = dims.ld;
155+ let xk = mat[k * ld + dims.n] / mat[k * ld + k];
156+ if (lid.x == 0u) { sol[k] = xk; }
157+ for (var i = lid.x; i < k; i = i + ${LANES}u) {
158+ mat[i * ld + dims.n] = mat[i * ld + dims.n] - mat[i * ld + k] * xk;
159+ }
160+}
161+
162+/** The potential at the targets: one thread per target, summing the
163+ * charges. */
164+@compute @workgroup_size(${EVAL_GROUP})
165+fn evaluate(@builtin(global_invocation_id) gid: vec3<u32>) {
166+ let t = gid.x;
167+ if (t >= dims.m) { return; }
168+ let ax = tg_x(t);
169+ let ay = tg_y(t);
170+ var acc = 0.0;
171+ for (var j = 0u; j < dims.n; j = j + 1u) {
172+ acc = acc + sol[j] * logdist(ax, ay, ch_x(j), ch_y(j));
173+ }
174+ out[t] = acc;
175+}
176+`;
177+}
src/solvers/webgpuSolvers.tsadded+45−0View file
@@ -0,0 +1,45 @@
1+// The WebGPU solvers, by manifest id.
2+//
3+// A solver whose runtime is "webgpu" is not a MATLAB file, so it cannot be
4+// looked up on disk the way the numbl and MATLAB ones are; it is a module
5+// in this repository that implements the small interface below against the
6+// TypeScript form of the problem (src/problems/laplace2d/problem.ts). One
7+// instance is built per id and kept: the shader compilation and pipeline
8+// creation it does belong outside the timed runs.
9+
10+import type { Laplace2dProblem } from "../problems/laplace2d/problem";
11+import { MfsGpu } from "./mfs-gpu/solver";
12+
13+export interface WebgpuSolver {
14+ /** One full solve at resolution n, ending when the device has finished
15+ * and the values have been read back. */
16+ run(
17+ prob: Laplace2dProblem,
18+ n: number,
19+ wantGrid: boolean
20+ ): Promise<{ uEval: Float64Array; uGrid: Float64Array | null }>;
21+ /** The adapter and how WebGPU was reached, for the result file. */
22+ readonly adapter: string;
23+ readonly via: string;
24+}
25+
26+const factories: Record<string, () => Promise<WebgpuSolver>> = {
27+ "mfs-gpu": () => MfsGpu.create(),
28+};
29+
30+const built = new Map<string, Promise<WebgpuSolver>>();
31+
32+export function isWebgpuSolver(id: string): boolean {
33+ return id in factories;
34+}
35+
36+export function getWebgpuSolver(id: string): Promise<WebgpuSolver> {
37+ const make = factories[id];
38+ if (!make) throw new Error(`no WebGPU solver named ${id}`);
39+ let p = built.get(id);
40+ if (!p) {
41+ p = make();
42+ built.set(id, p);
43+ }
44+ return p;
45+}
test/expected.tsmodified+17−0View file
@@ -24,6 +24,18 @@ export const MUST_REACH: Record<string, Record<string, number>> = {
2424 // nothing extra because the representation is smooth up to it.
2525 "star-nearfield": 1e-4,
2626 },
27+ // WebGPU has no f64, and the MFS is conditioning-limited, so these
28+ // floors are six to nine orders looser than the same method's on the CPU.
29+ // That is the finding, not a defect; MUST_NOT_REACH below holds it in
30+ // place.
31+ "mfs-gpu": {
32+ "disk-easy": 1e-6,
33+ "star-medium": 1e-6,
34+ "star-hard": 1e-2,
35+ "flower-15": 1e-2,
36+ "square-corners": 1e-5,
37+ "star-nearfield": 1e-2,
38+ },
2739 "nystrom-dlp": {
2840 "disk-easy": 1e-10,
2941 "star-medium": 1e-10,
@@ -65,4 +77,9 @@ export const MUST_NOT_REACH: Record<string, Record<string, number>> = {
6577 // a target 0.005 inside the boundary, either it acquired a near-field
6678 // correction or the instance stopped placing its targets there.
6779 "nystrom-dlp": { "star-nearfield": 1e-8 },
80+ // Single precision is what caps mfs-gpu on the two instances where the
81+ // method itself would otherwise reach 1e-13. If it ever got past this,
82+ // it stopped computing in f32 and the pair with mfs stopped measuring
83+ // what it claims to measure.
84+ "mfs-gpu": { "disk-easy": 1e-9, "star-medium": 1e-9 },
6885 };
test/gpu-test.tsadded+76−0View file
@@ -0,0 +1,76 @@
1+// Convergence test for the WebGPU-runtime solvers, run where a WebGPU
2+// device can be had: npx tsx test/gpu-test.ts
3+// Exits quietly with a notice when there is none, which is the usual case
4+// on a headless machine and in CI.
5+//
6+// The expectations live in test/expected.ts, keyed by source directory, as
7+// for the other two suites. mfs-gpu's floors are deliberately far looser
8+// than mfs's, because WebGPU has no f64; MUST_NOT_REACH pins that down from
9+// the other side, so a run that suddenly got double precision fails here
10+// rather than quietly changing what the pair of curves means.
11+
12+import { getInstance, INSTANCES } from "../src/problems/laplace2d/spec";
13+import { SOLVERS, solverSourceDir, sweepNFor } from "../src/solvers";
14+import { runSweepGpu } from "../src/harness/webgpuRun";
15+import { gpuUnavailableReason, requestGpu } from "../src/harness/webgpuDevice";
16+import { DEFAULT_TIMING } from "../src/harness/timing";
17+import { MUST_NOT_REACH, MUST_REACH } from "./expected";
18+
19+const why = await gpuUnavailableReason();
20+if (why !== null) {
21+ console.log(`no WebGPU device here; skipping the WebGPU solver tests\n ${why}`);
22+ process.exit(0);
23+}
24+const { adapter, via } = await requestGpu();
25+console.log(`WebGPU: ${adapter} (${via})`);
26+
27+// Accuracy, not speed: one timed run per point is enough here.
28+const TEST_TIMING = { ...DEFAULT_TIMING, minTimedRuns: 1, timeBudgetSeconds: 0 };
29+
30+let failures = 0;
31+for (const solver of SOLVERS.filter((s) => s.runtime === "webgpu")) {
32+ const dir = solverSourceDir(solver);
33+ for (const inst of INSTANCES) {
34+ const reach = MUST_REACH[dir]?.[inst.id];
35+ if (reach === undefined) continue;
36+ console.log(`\n== ${inst.id} / ${solver.id} (WebGPU)`);
37+ console.log(" n relMax relL2 solve(s)");
38+ let best = Infinity;
39+ const points = await runSweepGpu({
40+ instance: getInstance(inst.id),
41+ solver,
42+ timing: TEST_TIMING,
43+ onPoint: (p) => {
44+ best = Math.min(best, p.relMax);
45+ console.log(
46+ ` ${String(p.n).padStart(4)} ${p.relMax.toExponential(3)} ` +
47+ `${p.relL2.toExponential(3)} ${p.solveSeconds.toFixed(4)}`
48+ );
49+ },
50+ });
51+ if (points.length !== sweepNFor(solver, inst.id).length) {
52+ console.log(" FAIL: sweep returned the wrong number of points");
53+ failures++;
54+ continue;
55+ }
56+ const notReach = MUST_NOT_REACH[dir]?.[inst.id];
57+ if (best > reach) {
58+ console.log(` FAIL: best relMax ${best.toExponential(2)} > ${reach}`);
59+ failures++;
60+ } else if (notReach !== undefined && best < notReach) {
61+ console.log(
62+ ` FAIL: best relMax ${best.toExponential(2)} < ${notReach} ` +
63+ "(single precision no longer caps this solver)"
64+ );
65+ failures++;
66+ } else {
67+ console.log(` ok (best relMax ${best.toExponential(2)})`);
68+ }
69+ }
70+}
71+
72+if (failures > 0) {
73+ console.error(`\n${failures} failure(s)`);
74+ process.exit(1);
75+}
76+console.log("\nall WebGPU checks passed");
test/solver-test.tsmodified+13−8View file
@@ -12,6 +12,7 @@ import { INSTANCES, getInstance } from "../src/problems/laplace2d/spec";
1212 import {
1313 SOLVERS,
1414 getSolver,
15+ solverFiles,
1516 solverSourceDir,
1617 type SolverManifest,
1718 } from "../src/solvers";
@@ -45,22 +46,26 @@ const TEST_TIMING = { ...DEFAULT_TIMING, minTimedRuns: 1, timeBudgetSeconds: 0 }
4546 let failures = 0;
4647
4748 // Registry consistency, checked for every entry including the ones this
48-// suite does not run: the solver file must exist, it must have stated
49-// expectations, and an entry that borrows another entry's solver.m must
50-// carry the same version, so that two results claiming the same solver
49+// suite does not run: the solver's source files must exist, it must have
50+// stated expectations, and an entry that borrows another entry's source
51+// must carry the same version, so that two results claiming the same solver
5152 // version really did run the same code.
5253 for (const s of SOLVERS) {
5354 const dir = solverSourceDir(s);
54- const path = `src/solvers/${dir}/solver.m`;
55- if (!existsSync(join(root, path))) {
56- console.log(`FAIL: ${s.id} has no ${path}`);
57- failures++;
55+ for (const file of solverFiles(s)) {
56+ const path = `src/solvers/${file}`;
57+ if (!existsSync(join(root, path))) {
58+ console.log(`FAIL: ${s.id} has no ${path}`);
59+ failures++;
60+ }
5861 }
5962 const expected = MUST_REACH[dir];
6063 if (!expected) {
6164 console.log(`FAIL: ${s.id} has no entry in test/expected.ts`);
6265 failures++;
63- } else if (s.runtime === "numbl") {
66+ } else if (s.runtime === "numbl" || s.runtime === "webgpu") {
67+ // Both of these run every instance, in test/solver-test.ts and
68+ // test/gpu-test.ts respectively.
6469 for (const inst of INSTANCES) {
6570 if (expected[inst.id] === undefined) {
6671 console.log(`FAIL: ${s.id} has no expectation on ${inst.id}`);
tsconfig.jsonmodified+1−1View file
@@ -11,7 +11,7 @@
1111 "isolatedModules": true,
1212 "noUnusedLocals": true,
1313 "noUnusedParameters": true,
14- "types": ["vite/client", "node"]
14+ "types": ["vite/client", "node", "@webgpu/types"]
1515 },
1616 "include": ["src", "test", "scripts", "vite.config.ts"]
1717 }