fastandaccurate: PDE solver benchmarks, laplace-dirichlet-2d with MFS and Nystrom DLP solvers
41 changed files+7142−0
.github/workflows/deploy.ymladded+37−0View file
@@ -0,0 +1,37 @@
1+name: deploy
2+on:
3+ push:
4+ branches: [main]
5+ workflow_dispatch:
6+
7+permissions:
8+ contents: read
9+ pages: write
10+ id-token: write
11+
12+concurrency:
13+ group: pages
14+ cancel-in-progress: true
15+
16+jobs:
17+ build-deploy:
18+ runs-on: ubuntu-latest
19+ environment:
20+ name: github-pages
21+ url: ${{ steps.deployment.outputs.page_url }}
22+ steps:
23+ - uses: actions/checkout@v4
24+ - uses: actions/setup-node@v4
25+ with:
26+ node-version: 24
27+ cache: npm
28+ - run: npm ci
29+ - run: npm test
30+ - run: npm run build
31+ # Pages must already be enabled with "GitHub Actions" as the source.
32+ - uses: actions/configure-pages@v5
33+ - uses: actions/upload-pages-artifact@v3
34+ with:
35+ path: dist
36+ - id: deployment
37+ uses: actions/deploy-pages@v4
.gitignoreadded+5−0View file
@@ -0,0 +1,5 @@
1+node_modules/
2+dist/
3+dist-cli/
4+fastandaccurate-results-out/
5+*.tsbuildinfo
LICENSEadded+202−0View file
@@ -0,0 +1,202 @@
1+
2+ Apache License
3+ Version 2.0, January 2004
4+ http://www.apache.org/licenses/
5+
6+ TERMS AND CONDITIONS FOR USE, REPRODUCTION, AND DISTRIBUTION
7+
8+ 1. Definitions.
9+
10+ "License" shall mean the terms and conditions for use, reproduction,
11+ and distribution as defined by Sections 1 through 9 of this document.
12+
13+ "Licensor" shall mean the copyright owner or entity authorized by
14+ the copyright owner that is granting the License.
15+
16+ "Legal Entity" shall mean the union of the acting entity and all
17+ other entities that control, are controlled by, or are under common
18+ control with that entity. For the purposes of this definition,
19+ "control" means (i) the power, direct or indirect, to cause the
20+ direction or management of such entity, whether by contract or
21+ otherwise, or (ii) ownership of fifty percent (50%) or more of the
22+ outstanding shares, or (iii) beneficial ownership of such entity.
23+
24+ "You" (or "Your") shall mean an individual or Legal Entity
25+ exercising permissions granted by this License.
26+
27+ "Source" form shall mean the preferred form for making modifications,
28+ including but not limited to software source code, documentation
29+ source, and configuration files.
30+
31+ "Object" form shall mean any form resulting from mechanical
32+ transformation or translation of a Source form, including but
33+ not limited to compiled object code, generated documentation,
34+ and conversions to other media types.
35+
36+ "Work" shall mean the work of authorship, whether in Source or
37+ Object form, made available under the License, as indicated by a
38+ copyright notice that is included in or attached to the work
39+ (an example is provided in the Appendix below).
40+
41+ "Derivative Works" shall mean any work, whether in Source or Object
42+ form, that is based on (or derived from) the Work and for which the
43+ editorial revisions, annotations, elaborations, or other modifications
44+ represent, as a whole, an original work of authorship. For the purposes
45+ of this License, Derivative Works shall not include works that remain
46+ separable from, or merely link (or bind by name) to the interfaces of,
47+ the Work and Derivative Works thereof.
48+
49+ "Contribution" shall mean any work of authorship, including
50+ the original version of the Work and any modifications or additions
51+ to that Work or Derivative Works thereof, that is intentionally
52+ submitted to Licensor for inclusion in the Work by the copyright owner
53+ or by an individual or Legal Entity authorized to submit on behalf of
54+ the copyright owner. For the purposes of this definition, "submitted"
55+ means any form of electronic, verbal, or written communication sent
56+ to the Licensor or its representatives, including but not limited to
57+ communication on electronic mailing lists, source code control systems,
58+ and issue tracking systems that are managed by, or on behalf of, the
59+ Licensor for the purpose of discussing and improving the Work, but
60+ excluding communication that is conspicuously marked or otherwise
61+ designated in writing by the copyright owner as "Not a Contribution."
62+
63+ "Contributor" shall mean Licensor and any individual or Legal Entity
64+ on behalf of whom a Contribution has been received by Licensor and
65+ subsequently incorporated within the Work.
66+
67+ 2. Grant of Copyright License. Subject to the terms and conditions of
68+ this License, each Contributor hereby grants to You a perpetual,
69+ worldwide, non-exclusive, no-charge, royalty-free, irrevocable
70+ copyright license to reproduce, prepare Derivative Works of,
71+ publicly display, publicly perform, sublicense, and distribute the
72+ Work and such Derivative Works in Source or Object form.
73+
74+ 3. Grant of Patent License. Subject to the terms and conditions of
75+ this License, each Contributor hereby grants to You a perpetual,
76+ worldwide, non-exclusive, no-charge, royalty-free, irrevocable
77+ (except as stated in this section) patent license to make, have made,
78+ use, offer to sell, sell, import, and otherwise transfer the Work,
79+ where such license applies only to those patent claims licensable
80+ by such Contributor that are necessarily infringed by their
81+ Contribution(s) alone or by combination of their Contribution(s)
82+ with the Work to which such Contribution(s) was submitted. If You
83+ institute patent litigation against any entity (including a
84+ cross-claim or counterclaim in a lawsuit) alleging that the Work
85+ or a Contribution incorporated within the Work constitutes direct
86+ or contributory patent infringement, then any patent licenses
87+ granted to You under this License for that Work shall terminate
88+ as of the date such litigation is filed.
89+
90+ 4. Redistribution. You may reproduce and distribute copies of the
91+ Work or Derivative Works thereof in any medium, with or without
92+ modifications, and in Source or Object form, provided that You
93+ meet the following conditions:
94+
95+ (a) You must give any other recipients of the Work or
96+ Derivative Works a copy of this License; and
97+
98+ (b) You must cause any modified files to carry prominent notices
99+ stating that You changed the files; and
100+
101+ (c) You must retain, in the Source form of any Derivative Works
102+ that You distribute, all copyright, patent, trademark, and
103+ attribution notices from the Source form of the Work,
104+ excluding those notices that do not pertain to any part of
105+ the Derivative Works; and
106+
107+ (d) If the Work includes a "NOTICE" text file as part of its
108+ distribution, then any Derivative Works that You distribute must
109+ include a readable copy of the attribution notices contained
110+ within such NOTICE file, excluding those notices that do not
111+ pertain to any part of the Derivative Works, in at least one
112+ of the following places: within a NOTICE text file distributed
113+ as part of the Derivative Works; within the Source form or
114+ documentation, if provided along with the Derivative Works; or,
115+ within a display generated by the Derivative Works, if and
116+ wherever such third-party notices normally appear. The contents
117+ of the NOTICE file are for informational purposes only and
118+ do not modify the License. You may add Your own attribution
119+ notices within Derivative Works that You distribute, alongside
120+ or as an addendum to the NOTICE text from the Work, provided
121+ that such additional attribution notices cannot be construed
122+ as modifying the License.
123+
124+ You may add Your own copyright statement to Your modifications and
125+ may provide additional or different license terms and conditions
126+ for use, reproduction, or distribution of Your modifications, or
127+ for any such Derivative Works as a whole, provided Your use,
128+ reproduction, and distribution of the Work otherwise complies with
129+ the conditions stated in this License.
130+
131+ 5. Submission of Contributions. Unless You explicitly state otherwise,
132+ any Contribution intentionally submitted for inclusion in the Work
133+ by You to the Licensor shall be under the terms and conditions of
134+ this License, without any additional terms or conditions.
135+ Notwithstanding the above, nothing herein shall supersede or modify
136+ the terms of any separate license agreement you may have executed
137+ with Licensor regarding such Contributions.
138+
139+ 6. Trademarks. This License does not grant permission to use the trade
140+ names, trademarks, service marks, or product names of the Licensor,
141+ except as required for reasonable and customary use in describing the
142+ origin of the Work and reproducing the content of the NOTICE file.
143+
144+ 7. Disclaimer of Warranty. Unless required by applicable law or
145+ agreed to in writing, Licensor provides the Work (and each
146+ Contributor provides its Contributions) on an "AS IS" BASIS,
147+ WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or
148+ implied, including, without limitation, any warranties or conditions
149+ of TITLE, NON-INFRINGEMENT, MERCHANTABILITY, or FITNESS FOR A
150+ PARTICULAR PURPOSE. You are solely responsible for determining the
151+ appropriateness of using or redistributing the Work and assume any
152+ risks associated with Your exercise of permissions under this License.
153+
154+ 8. Limitation of Liability. In no event and under no legal theory,
155+ whether in tort (including negligence), contract, or otherwise,
156+ unless required by applicable law (such as deliberate and grossly
157+ negligent acts) or agreed to in writing, shall any Contributor be
158+ liable to You for damages, including any direct, indirect, special,
159+ incidental, or consequential damages of any character arising as a
160+ result of this License or out of the use or inability to use the
161+ Work (including but not limited to damages for loss of goodwill,
162+ work stoppage, computer failure or malfunction, or any and all
163+ other commercial damages or losses), even if such Contributor
164+ has been advised of the possibility of such damages.
165+
166+ 9. Accepting Warranty or Additional Liability. While redistributing
167+ the Work or Derivative Works thereof, You may choose to offer,
168+ and charge a fee for, acceptance of support, warranty, indemnity,
169+ or other liability obligations and/or rights consistent with this
170+ License. However, in accepting such obligations, You may act only
171+ on Your own behalf and on Your sole responsibility, not on behalf
172+ of any other Contributor, and only if You agree to indemnify,
173+ defend, and hold each Contributor harmless for any liability
174+ incurred by, or claims asserted against, such Contributor by reason
175+ of your accepting any such warranty or additional liability.
176+
177+ END OF TERMS AND CONDITIONS
178+
179+ APPENDIX: How to apply the Apache License to your work.
180+
181+ To apply the Apache License to your work, attach the following
182+ boilerplate notice, with the fields enclosed by brackets "[]"
183+ replaced with your own identifying information. (Don't include
184+ the brackets!) The text should be enclosed in the appropriate
185+ comment syntax for the file format. We also recommend that a
186+ file or class name and description of purpose be included on the
187+ same "printed page" as the copyright notice for easier
188+ identification within third-party archives.
189+
190+ Copyright [yyyy] [name of copyright owner]
191+
192+ Licensed under the Apache License, Version 2.0 (the "License");
193+ you may not use this file except in compliance with the License.
194+ You may obtain a copy of the License at
195+
196+ http://www.apache.org/licenses/LICENSE-2.0
197+
198+ Unless required by applicable law or agreed to in writing, software
199+ distributed under the License is distributed on an "AS IS" BASIS,
200+ WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
201+ See the License for the specific language governing permissions and
202+ limitations under the License.
README.mdadded+87−0View file
@@ -0,0 +1,87 @@
1+# fastandaccurate
2+
3+Speed and accuracy benchmarks for PDE solvers.
4+
5+Live site: https://concept-collection.github.io/fastandaccurate/
6+
7+A limitation of most solver comparisons is that they fix a
8+discretization, which quietly decides much of the outcome. Here each
9+**problem** is posed in the continuum with an exact reference solution;
10+a solver chooses its own discretization and is scored at
11+problem-specified evaluation points. The central object is the
12+**work-precision curve**: error against compute time as the solver's
13+resolution varies. There is deliberately no single ranking, since which
14+curve wins can differ by accuracy regime, instance, and machine.
15+
16+Solvers are MATLAB function files run by [numbl](https://numbl.org)
17+(MATLAB syntax in the browser and in node), so every run on this site
18+happens client side, and the identical harness runs from the command
19+line. Each problem defines its own interface and instances in a written
20+specification; the system is deliberately loose, and per-problem
21+interfaces are expected to differ.
22+
23+## Problems
24+
25+- [laplace-dirichlet-2d](docs/problems/laplace-dirichlet-2d.md) —
26+ interior Dirichlet Laplace problem on a star-shaped domain, data
27+ manufactured from an exact harmonic function whose singularities sit an
28+ adjustable distance outside the boundary.
29+
30+## Results
31+
32+Results are work-precision sweeps stored as JSON files in
33+[fastandaccurate-results](https://github.com/concept-collection/fastandaccurate-results)
34+and added by pull request; the site reads that repository statically, so
35+there is no database and no server. Every result records its provenance:
36+instance spec and hash, solver id and version, protocol, runtime, numbl
37+version, and machine. In-browser results can be rerun by any visitor on
38+their own machine directly on the site.
39+
40+## Running benchmarks outside the browser
41+
42+The command line installs from the site itself (node 20 or newer):
43+
44+```
45+npx https://concept-collection.github.io/fastandaccurate/cli.tgz run --label "my workstation"
46+```
47+
48+Note that npx caches by the exact URL string; the site offers the URL
49+with a `?v=<commit>` suffix so each deployment is a fresh install.
50+
51+Useful flags: `--instance <id>`, `--solver <id>`, `--repeats N`,
52+`--max-n N`, `--out dir`. To benchmark your own solver, point the
53+harness at a MATLAB function file implementing the problem's interface:
54+
55+```
56+npx https://concept-collection.github.io/fastandaccurate/cli.tgz run \
57+ --solver-file my_method.m --solver-id my-method
58+```
59+
60+The resulting JSON files can be loaded on the site (load result file) to
61+view them against the committed curves, and submitted by PR to the
62+results repository. To add a solver to the site itself (so visitors can
63+rerun it in the browser), PR the solver directory and a manifest entry
64+to this repository; see `src/solvers/`.
65+
66+## Development
67+
68+```
69+npm install
70+npm run dev # local dev server
71+npm test # solver convergence tests through numbl in node
72+npm run build # type-check, site build, CLI tarball (dist/)
73+npm run check-app # headless end-to-end check of the built site
74+```
75+
76+Layout: `src/problems/` holds problem specs, instances, exact solutions,
77+and the problem-side MATLAB; `src/solvers/` the solver MATLAB files and
78+manifests; `src/harness/` the shared runner, sweep, and result schema
79+(used identically by the browser worker and the CLI); `src/app/` the
80+React site; `src/cli/` the command line.
81+
82+Deployed to GitHub Pages by `.github/workflows/deploy.yml` on push to
83+main.
84+
85+## License
86+
87+Apache-2.0
docs/problems/laplace-dirichlet-2d.mdadded+120−0View file
@@ -0,0 +1,120 @@
1+# Problem: laplace-dirichlet-2d (version 1)
2+
3+This document is the canonical statement of the problem. Results refer to
4+it by the pair (problem id, problem version); any change that could alter
5+a measured number requires a version bump.
6+
7+## The problem
8+
9+Solve the Laplace equation
10+
11+ Δu = 0 in Ω, u = g on ∂Ω,
12+
13+where Ω is the star-shaped plane domain bounded by the curve
14+
15+ x(t) = r(t) (cos t, sin t), r(t) = 1 + a cos(k t), t ∈ [0, 2π),
16+
17+and g is Dirichlet data specified below. The problem is posed in the
18+continuum: a solver receives the curve and the data as functions, chooses
19+its own discretization, and is scored on the values it returns at a fixed
20+set of evaluation points.
21+
22+## The exact solution
23+
24+The data g is manufactured from an exact harmonic function, a sum of
25+three logarithmic point sources placed outside the domain:
26+
27+ u*(x) = Σ_j c_j log |x − s_j|, c = (1.0, −0.6, 0.8).
28+
29+Source s_j is the boundary point at parameter φ_j = 2π(j−1)/3 + 0.4
30+pushed a distance d along the outward unit normal, and g = u* restricted
31+to ∂Ω. Since u* is harmonic in Ω, it is the unique solution, and errors
32+are measured against it directly rather than against a reference
33+computation. The distance d controls difficulty: u* continues
34+harmonically only up to the sources, so the smaller d, the shorter the
35+distance the data continues past the boundary, and methods whose
36+representations assume a generous continuation lose it. A solver must
37+not use knowledge of the sources; they exist only to manufacture g.
38+
39+## Official instances
40+
41+| id | a | k | d | character |
42+|---|---|---|---|---|
43+| disk-easy | 0 | 0 | 0.5 | the unit disk, distant sources |
44+| star-medium | 0.2 | 3 | 0.4 | mild geometry, comfortable continuation |
45+| star-hard | 0.3 | 5 | 0.08 | wavy boundary, data barely continues |
46+
47+Committed results exist only at these instances, so every solver is
48+compared on identical inputs. The parameters may be varied freely in the
49+site's interactive views, but such runs are not recorded.
50+
51+## What is scored
52+
53+Each instance fixes 65 evaluation points: 16 rays at angles
54+θ_j = 2πj/16 + 0.13 (j = 0, …, 15), radial fractions
55+ρ ∈ {0.25, 0.5, 0.75, 0.9} along each ray (the point at fraction ρ on ray
56+θ is ρ r(θ) (cos θ, sin θ)), ordered radius-major, plus the origin last.
57+The solver returns u at exactly these points. Reported errors are
58+
59+ relMax = max_i |u_i − u*_i| / max_i |u*_i|,
60+ relL2 = ( Σ_i (u_i − u*_i)² / Σ_i u*_i² )^{1/2}.
61+
62+relMax is the headline number; both are recorded. Note that the
63+evaluation points deliberately stay a modest distance inside the domain
64+(fraction 0.9 at most). Accuracy very close to the boundary is a genuine
65+and separate difficulty (the close-evaluation problem for integral
66+methods, boundary-layer resolution for grid methods) and a future
67+problem version may add a near-boundary target set; version 1 does not
68+test it.
69+
70+## Timing protocol
71+
72+One run is a full call of the solver, including its own discretization,
73+assembly, solve, and evaluation. The harness performs one untimed warmup
74+run (which absorbs numbl's JIT compilation and is recorded separately as
75+the cold time), then N timed runs (N = 3 unless stated), and reports the
76+median as the solve time. All timing is tic/toc inside the MATLAB
77+session, so browser and node runs measure the same thing. Times from
78+different machines are not comparable; every result records its
79+environment, and comparisons across environments are the reader's
80+responsibility.
81+
82+## Solver interface
83+
84+A solver is a MATLAB function file
85+
86+ function out = solver(prob, n)
87+
88+where n is the solver's own resolution parameter (its meaning is the
89+solver's choice; the standard sweep list is declared in the solver's
90+manifest) and prob is a struct with fields
91+
92+| field | meaning |
93+|---|---|
94+| `prob.curve` | `@(t) -> [x y]`, boundary point at parameter t (column vectors in, m×2 out) |
95+| `prob.curveD` | first derivative of the curve with respect to t |
96+| `prob.curveDD` | second derivative |
97+| `prob.g` | `@(t) -> g`, Dirichlet data at boundary parameter t |
98+| `prob.evalXY` | 65×2, the evaluation points |
99+| `prob.vizXY` | m×2, visualization grid points (m = 0 when not requested) |
100+
101+The return value is a struct: `out.uEval` (65×1, required) and
102+`out.uGrid` (m×1; `[]` when `prob.vizXY` is empty). The solver must not
103+reconstruct the sources analytically or otherwise special-case the known
104+solution; submissions are reviewed for this.
105+
106+## Visualization grid
107+
108+When requested, `prob.vizXY` lists a 200×200 grid of points over the
109+bounding square [−R, R]², R = 1.05(1 + |a|), with flat index
110+p = ix·200 + iy for x = xs[ix], y = xs[iy] (y varies fastest, MATLAB
111+meshgrid column order). Points outside Ω are included and the viewer
112+masks them; grid values are never scored.
113+
114+## Limitations
115+
116+The exact solution is smooth and free of boundary singularities, so this
117+problem does not test corner handling, nonsmooth data, or interior
118+sources (a nonzero right-hand side would exclude plain boundary-integral
119+methods; that belongs to a different problem). The domain family is
120+star-shaped by construction, which some methods can exploit.
index.htmladded+20−0View file
@@ -0,0 +1,20 @@
1+<!doctype html>
2+<html lang="en">
3+ <head>
4+ <meta charset="UTF-8" />
5+ <meta name="viewport" content="width=device-width, initial-scale=1.0" />
6+ <title>fastandaccurate — PDE solver benchmarks</title>
7+ <link
8+ rel="icon"
9+ href="data:image/svg+xml,%3Csvg xmlns='http://www.w3.org/2000/svg' viewBox='0 0 16 16'%3E%3Ctext y='13' font-size='13'%3E%E2%88%86%3C/text%3E%3C/svg%3E"
10+ />
11+ <meta
12+ name="description"
13+ content="Speed and accuracy benchmarks for PDE solvers: continuum problems with exact ground truth, work-precision comparisons, reproducible in the browser."
14+ />
15+ </head>
16+ <body>
17+ <div id="root"></div>
18+ <script type="module" src="/src/app/main.tsx"></script>
19+ </body>
20+</html>
package-lock.jsonadded+3626−0View file
This diff is 3,631 lines long and is not shown.
package.jsonadded+28−0View file
@@ -0,0 +1,28 @@
1+{
2+ "name": "fastandaccurate",
3+ "private": true,
4+ "version": "0.1.0",
5+ "type": "module",
6+ "scripts": {
7+ "dev": "vite",
8+ "build": "tsc --noEmit && vite build && npm run build:cli",
9+ "build:cli": "vite build --config vite.cli.config.ts && node scripts/pack-cli.mjs",
10+ "preview": "vite preview",
11+ "test": "tsx test/solver-test.ts",
12+ "check-app": "node scripts/check-app.mjs"
13+ },
14+ "dependencies": {
15+ "numbl": "^0.4.18",
16+ "react": "^19.2.7",
17+ "react-dom": "^19.2.7"
18+ },
19+ "devDependencies": {
20+ "@types/node": "^24.0.0",
21+ "@types/react": "^19.2.0",
22+ "@types/react-dom": "^19.2.0",
23+ "@vitejs/plugin-react": "^5.0.0",
24+ "tsx": "^4.19.0",
25+ "typescript": "~5.9.3",
26+ "vite": "^7.0.0"
27+ }
28+}
scripts/check-app.mjsadded+121−0View file
@@ -0,0 +1,121 @@
1+/**
2+ * Headless end-to-end check of the built site: serves dist/, loads the
3+ * page in headless Chrome, fails on any console error, and exercises one
4+ * real in-browser solve through the worker (the Solution section),
5+ * asserting the reported error is small. Visual appearance is checked by
6+ * a human, not here.
7+ *
8+ * Usage: npm run build && node scripts/check-app.mjs
9+ */
10+import { createServer } from "node:http";
11+import { readFile } from "node:fs/promises";
12+import { join, extname } from "node:path";
13+import puppeteer from "puppeteer-core";
14+
15+const root = new URL("../dist", import.meta.url).pathname;
16+const types = {
17+ ".html": "text/html",
18+ ".js": "text/javascript",
19+ ".css": "text/css",
20+ ".wasm": "application/wasm",
21+ ".json": "application/json",
22+};
23+
24+const server = createServer(async (req, res) => {
25+ const path = req.url === "/" ? "/index.html" : (req.url ?? "/").split("?")[0];
26+ try {
27+ const data = await readFile(join(root, path));
28+ res.writeHead(200, { "content-type": types[extname(path)] ?? "application/octet-stream" });
29+ res.end(data);
30+ } catch {
31+ res.writeHead(404);
32+ res.end("not found");
33+ }
34+});
35+await new Promise((resolve) => server.listen(0, resolve));
36+const port = server.address().port;
37+
38+const browser = await puppeteer.launch({
39+ executablePath: process.env.CHROME_PATH ?? "/usr/bin/google-chrome",
40+ args: ["--no-sandbox"],
41+});
42+
43+const errors = [];
44+let failures = 0;
45+try {
46+ const page = await browser.newPage();
47+ page.on("console", (msg) => {
48+ if (msg.type() === "error") {
49+ const text = msg.text();
50+ const url = msg.location()?.url ?? "";
51+ // The committed-results fetch may 404 before the results repo has
52+ // content; the app reports that in the UI by design.
53+ if (url.includes("raw.githubusercontent.com")) return;
54+ // numbl logs its linear-algebra bridge choice at error level.
55+ if (text.includes("using bridge:")) return;
56+ errors.push(`${text} (${url})`);
57+ }
58+ });
59+ page.on("pageerror", (err) => errors.push(String(err)));
60+
61+ await page.goto(`http://127.0.0.1:${port}/`, { waitUntil: "networkidle2" });
62+ await page.waitForSelector("h1");
63+
64+ // The empty-state or chart should be present.
65+ const hasChartArea = await page.evaluate(
66+ () => document.body.innerText.includes("Work-precision results")
67+ );
68+ if (!hasChartArea) {
69+ console.error("FAIL: work-precision section missing");
70+ failures++;
71+ }
72+
73+ // One real solve through the worker: the Solution section's compute
74+ // button, default solver (mfs) at its default n.
75+ const buttons = await page.$$("button");
76+ let computeBtn = null;
77+ for (const b of buttons) {
78+ const t = await b.evaluate((el) => el.textContent);
79+ if (t && t.includes("Compute in this browser")) computeBtn = b;
80+ }
81+ if (!computeBtn) {
82+ console.error("FAIL: compute button not found");
83+ failures++;
84+ } else {
85+ await computeBtn.click();
86+ await page.waitForFunction(
87+ () => document.body.innerText.includes("rel max error"),
88+ { timeout: 120000 }
89+ );
90+ const text = await page.evaluate(() => document.body.innerText);
91+ const m = text.match(/rel max error ([0-9.]+e[+-][0-9]+)/);
92+ if (!m) {
93+ console.error("FAIL: no reported error after compute");
94+ failures++;
95+ } else {
96+ const err = parseFloat(m[1]);
97+ // mfs at its default n on star-medium should be far below 1e-6.
98+ if (!(err < 1e-6)) {
99+ console.error(`FAIL: in-browser mfs error ${err} not < 1e-6`);
100+ failures++;
101+ } else {
102+ console.log(`in-browser solve ok (rel max error ${m[1]})`);
103+ }
104+ }
105+ }
106+
107+ if (errors.length > 0) {
108+ console.error("FAIL: console errors:");
109+ for (const e of errors) console.error(` ${e}`);
110+ failures++;
111+ }
112+} finally {
113+ await browser.close();
114+ server.close();
115+}
116+
117+if (failures > 0) {
118+ console.error(`${failures} failure(s)`);
119+ process.exit(1);
120+}
121+console.log("check-app: all checks passed");
scripts/cli-launcher.cjsadded+33−0View file
@@ -0,0 +1,33 @@
1+#!/usr/bin/env node
2+/**
3+ * The published command's entry point, whose only job is to be parseable
4+ * by whatever node it lands on. main.js is a modern bundle: on an old node
5+ * its first `import {` is a syntax error pointing at a brace, before any
6+ * check inside it could run. This file is therefore deliberately ES5, and
7+ * the dynamic import is built at run time, which older node cannot parse
8+ * as written syntax either.
9+ */
10+"use strict";
11+
12+var MIN_MAJOR = 20;
13+var have = process.versions.node;
14+var major = parseInt(have.split(".")[0], 10);
15+
16+if (!(major >= MIN_MAJOR)) {
17+ process.stderr.write(
18+ "fastandaccurate: this is node " + have + ", and the harness needs node " +
19+ MIN_MAJOR + " or newer (WebCrypto as a global, among other things).\n" +
20+ " nodejs.org has current builds; nvm, fnm and asdf install one per user\n" +
21+ " without touching what the system depends on.\n"
22+ );
23+ process.exit(1);
24+}
25+
26+var path = require("path");
27+var url = require("url");
28+var target = url.pathToFileURL(path.join(__dirname, "main.js")).href;
29+
30+new Function("specifier", "return import(specifier);")(target).catch(function (e) {
31+ process.stderr.write(String((e && e.stack) || e) + "\n");
32+ process.exit(1);
33+});
scripts/pack-cli.mjsadded+67−0View file
@@ -0,0 +1,67 @@
1+/**
2+ * Pack the command-line bundle as an npm tarball and put it in dist/,
3+ * where the page is deployed, so that
4+ *
5+ * npx https://concept-collection.github.io/fastandaccurate/cli.tgz
6+ *
7+ * installs and runs it anywhere node 20+ is. Nothing goes to the npm
8+ * registry: npm installs a tarball from a URL as happily as from a name,
9+ * and this way the command line is always the same commit as the site.
10+ * (npx caches by the exact URL string, so the site offers the URL with a
11+ * ?v=<commit> suffix to make each deployment a fresh install.)
12+ *
13+ * Usage: node scripts/pack-cli.mjs (after vite build --config vite.cli.config.ts)
14+ */
15+import { execFileSync } from "node:child_process";
16+import {
17+ mkdtemp,
18+ copyFile,
19+ writeFile,
20+ readFile,
21+ rm,
22+ mkdir,
23+} from "node:fs/promises";
24+import { glob } from "node:fs/promises";
25+import { tmpdir } from "node:os";
26+import { join, dirname, relative } from "node:path";
27+
28+const root = new URL("..", import.meta.url).pathname;
29+const pkg = JSON.parse(await readFile(join(root, "package.json"), "utf8"));
30+
31+const manifest = {
32+ name: "fastandaccurate-cli",
33+ version: pkg.version,
34+ description:
35+ "Run fastandaccurate PDE-solver benchmarks from the command line",
36+ license: "Apache-2.0",
37+ type: "module",
38+ engines: { node: ">=20" },
39+ bin: { fastandaccurate: "launch.cjs" },
40+};
41+
42+const stage = await mkdtemp(join(tmpdir(), "fastandaccurate-cli-"));
43+try {
44+ await copyFile(join(root, "dist-cli/main.js"), join(stage, "main.js"));
45+ await copyFile(join(root, "scripts/cli-launcher.cjs"), join(stage, "launch.cjs"));
46+ // The MATLAB sources, preserved under src/ so the CLI's source-root
47+ // probe finds them beside the bundle.
48+ for await (const entry of glob("src/{problems,solvers}/**/*.m", { cwd: root })) {
49+ const rel = relative(root, join(root, entry));
50+ await mkdir(join(stage, dirname(rel)), { recursive: true });
51+ await copyFile(join(root, rel), join(stage, rel));
52+ }
53+ await writeFile(
54+ join(stage, "package.json"),
55+ `${JSON.stringify(manifest, null, 2)}\n`
56+ );
57+ const out = execFileSync(
58+ "npm",
59+ ["pack", "--silent", "--pack-destination", stage],
60+ { cwd: stage, encoding: "utf8" }
61+ ).trim();
62+ await mkdir(join(root, "dist"), { recursive: true });
63+ await copyFile(join(stage, out), join(root, "dist/cli.tgz"));
64+ console.log(`dist/cli.tgz (${manifest.name} ${manifest.version})`);
65+} finally {
66+ await rm(stage, { recursive: true, force: true });
67+}
src/app/App.tsxadded+390−0View file
@@ -0,0 +1,390 @@
1+import { useEffect, useMemo, useState } from "react";
2+import {
3+ INSTANCES,
4+ getInstance,
5+ PROBLEM_ID,
6+} from "../problems/laplace2d/spec";
7+import { SOLVERS, getSolver } from "../solvers";
8+import {
9+ buildResultFile,
10+ type ResultFile,
11+ type ResultPoint,
12+} from "../harness/resultSchema";
13+import {
14+ environmentLabel,
15+ fetchCommittedResults,
16+ isResultFile,
17+ RESULTS_REPO_URL,
18+} from "./results";
19+import { solverColorVar } from "./colors";
20+import { sweepInBrowser } from "./workerClient";
21+import { WorkPrecisionChart, type ChartCurve } from "./components/WorkPrecisionChart";
22+import { PointsTable } from "./components/PointsTable";
23+import { DomainView } from "./components/DomainView";
24+import { SolutionSection } from "./components/SolutionSection";
25+
26+const REPO_URL = "https://github.com/concept-collection/fastandaccurate";
27+const SPEC_URL = `${REPO_URL}/blob/main/docs/problems/laplace-dirichlet-2d.md`;
28+
29+interface LocalRun {
30+ key: string;
31+ solverId: string;
32+ instanceId: string;
33+ repeats: number;
34+ points: ResultPoint[];
35+ done: boolean;
36+}
37+
38+export function App() {
39+ const [instanceId, setInstanceId] = useState(INSTANCES[1].id);
40+ const [committed, setCommitted] = useState<ResultFile[] | null>(null);
41+ const [committedError, setCommittedError] = useState<string | null>(null);
42+ const [loaded, setLoaded] = useState<ResultFile[]>([]);
43+ const [localRuns, setLocalRuns] = useState<LocalRun[]>([]);
44+ const [hidden, setHidden] = useState<Set<string>>(new Set());
45+ const [running, setRunning] = useState<string | null>(null);
46+ const [runStatus, setRunStatus] = useState<string | null>(null);
47+ const [repeats, setRepeats] = useState(3);
48+ const [machineLabel, setMachineLabel] = useState("");
49+
50+ const inst = getInstance(instanceId);
51+
52+ useEffect(() => {
53+ fetchCommittedResults()
54+ .then(setCommitted)
55+ .catch((err) =>
56+ setCommittedError(err instanceof Error ? err.message : String(err))
57+ );
58+ }, []);
59+
60+ const allSolverIds = useMemo(() => {
61+ const ids = new Set<string>(SOLVERS.map((s) => s.id));
62+ committed?.forEach((r) => ids.add(r.solver.id));
63+ loaded.forEach((r) => ids.add(r.solver.id));
64+ return [...ids];
65+ }, [committed, loaded]);
66+
67+ const curves: ChartCurve[] = useMemo(() => {
68+ const out: ChartCurve[] = [];
69+ const color = (id: string) => solverColorVar(id, allSolverIds);
70+ committed
71+ ?.filter((r) => r.problem === PROBLEM_ID && r.instance === instanceId)
72+ .forEach((r, i) => {
73+ out.push({
74+ key: `committed:${i}`,
75+ solverId: r.solver.id,
76+ label: `${r.solver.id} — ${environmentLabel(r)}`,
77+ color: color(r.solver.id),
78+ points: r.points,
79+ });
80+ });
81+ loaded
82+ .filter((r) => r.problem === PROBLEM_ID && r.instance === instanceId)
83+ .forEach((r, i) => {
84+ out.push({
85+ key: `loaded:${i}`,
86+ solverId: r.solver.id,
87+ label: `${r.solver.id} — ${environmentLabel(r)} (loaded)`,
88+ color: color(r.solver.id),
89+ dash: "8 4",
90+ points: r.points,
91+ });
92+ });
93+ localRuns
94+ .filter((r) => r.instanceId === instanceId)
95+ .forEach((r) => {
96+ out.push({
97+ key: r.key,
98+ solverId: r.solverId,
99+ label: `${r.solverId} — this browser`,
100+ color: color(r.solverId),
101+ dash: "4 4",
102+ open: true,
103+ points: r.points,
104+ });
105+ });
106+ return out.filter((c) => !hidden.has(c.solverId));
107+ }, [committed, loaded, localRuns, instanceId, hidden, allSolverIds]);
108+
109+ async function runSolver(solverId: string) {
110+ const key = `local:${solverId}:${instanceId}:${Date.now()}`;
111+ setLocalRuns((rs) => [
112+ ...rs.filter((r) => !(r.solverId === solverId && r.instanceId === instanceId)),
113+ { key, solverId, instanceId, repeats, points: [], done: false },
114+ ]);
115+ setRunning(solverId);
116+ try {
117+ const points = await sweepInBrowser(
118+ instanceId,
119+ solverId,
120+ repeats,
121+ (point, index, total) => {
122+ setRunStatus(
123+ `${solverId} on ${instanceId}: point ${index + 1}/${total} (n = ${point.n}) — rel max error ${point.relMax.toExponential(2)}`
124+ );
125+ setLocalRuns((rs) =>
126+ rs.map((r) => (r.key === key ? { ...r, points: [...r.points, point] } : r))
127+ );
128+ }
129+ );
130+ setLocalRuns((rs) =>
131+ rs.map((r) => (r.key === key ? { ...r, points, done: true } : r))
132+ );
133+ setRunStatus(null);
134+ } catch (err) {
135+ setRunStatus(
136+ `${solverId} failed: ${err instanceof Error ? err.message : String(err)}`
137+ );
138+ } finally {
139+ setRunning(null);
140+ }
141+ }
142+
143+ async function downloadRun(run: LocalRun) {
144+ const manifest = getSolver(run.solverId);
145+ const result = await buildResultFile({
146+ instance: getInstance(run.instanceId),
147+ solver: {
148+ id: manifest.id,
149+ version: manifest.version,
150+ backend: manifest.backend,
151+ source: "builtin",
152+ },
153+ environment: {
154+ kind: "browser",
155+ runtime: navigator.userAgent,
156+ numblVersion: __NUMBL_VERSION__,
157+ machineLabel: machineLabel || undefined,
158+ browserReproducible: true,
159+ },
160+ repeats: run.repeats,
161+ points: run.points,
162+ });
163+ const blob = new Blob([JSON.stringify(result, null, 2) + "\n"], {
164+ type: "application/json",
165+ });
166+ const a = document.createElement("a");
167+ a.href = URL.createObjectURL(blob);
168+ a.download = `${PROBLEM_ID}.${run.instanceId}.${run.solverId}.browser.json`;
169+ a.click();
170+ URL.revokeObjectURL(a.href);
171+ }
172+
173+ function loadFiles(files: FileList | null) {
174+ if (!files) return;
175+ for (const file of Array.from(files)) {
176+ file.text().then((text) => {
177+ try {
178+ const data: unknown = JSON.parse(text);
179+ if (isResultFile(data)) {
180+ setLoaded((ls) => [...ls, data]);
181+ } else {
182+ alert(`${file.name} is not a fastandaccurate result file`);
183+ }
184+ } catch {
185+ alert(`${file.name}: not valid JSON`);
186+ }
187+ });
188+ }
189+ }
190+
191+ const cliUrl = `https://concept-collection.github.io/fastandaccurate/cli.tgz?v=${__BUILD_ID__}`;
192+
193+ return (
194+ <main>
195+ <h1>fastandaccurate</h1>
196+ <p className="subtitle">Speed and accuracy benchmarks for PDE solvers</p>
197+ <p>
198+ Each <strong>problem</strong> here is posed in the continuum, with an
199+ exact reference solution; a solver chooses its own discretization and
200+ is scored at problem-specified evaluation points. The central object
201+ is the <strong>work-precision curve</strong>: error against compute
202+ time as the solver's resolution varies. There is deliberately no
203+ single ranking; which curve is best can differ by accuracy regime,
204+ instance, and machine. Results shown here are committed to a public{" "}
205+ <a href={RESULTS_REPO_URL}>results repository</a> by pull request, and
206+ any in-browser solver can be rerun on your own machine, right on this
207+ page, to check them.
208+ </p>
209+ <p className="small muted">
210+ <a href={REPO_URL}>Source</a> · <a href={SPEC_URL}>Problem specification</a> ·{" "}
211+ <a href={RESULTS_REPO_URL}>Results repository</a> · Solvers run in
212+ MATLAB syntax via <a href="https://numbl.org">numbl</a>, client side.
213+ </p>
214+
215+ <h2>The problem: laplace-dirichlet-2d</h2>
216+ <p>
217+ Solve Δu = 0 on the star-shaped domain with boundary
218+ r(θ) = 1 + a·cos(kθ), with Dirichlet data u = g on the boundary. The
219+ data comes from an exact harmonic function, a sum of three logarithmic
220+ point sources placed a distance d outside the boundary, so errors are
221+ measured against the true solution, not a reference computation. The
222+ distance d sets the difficulty: the closer the sources, the shorter
223+ the distance the data continues harmonically past the boundary, and
224+ methods whose representations assume that continuation lose it. A
225+ solver receives the curve (with derivatives), the boundary data as a
226+ function of the boundary parameter, and the evaluation points, and
227+ returns solution values at those points; reported time is the whole
228+ solve including the solver's own discretization (median of repeats
229+ after one untimed warmup). The precise statement, interface, and
230+ protocol are in the <a href={SPEC_URL}>specification</a>.
231+ </p>
232+ <div className="row" style={{ marginTop: 14 }}>
233+ <div>
234+ <div style={{ marginBottom: 10 }}>
235+ <label>
236+ instance{" "}
237+ <select
238+ value={instanceId}
239+ onChange={(e) => setInstanceId(e.target.value)}
240+ >
241+ {INSTANCES.map((i) => (
242+ <option key={i.id} value={i.id}>
243+ {i.id} — {i.label}
244+ </option>
245+ ))}
246+ </select>
247+ </label>
248+ </div>
249+ <p className="small muted" style={{ maxWidth: 380 }}>
250+ {inst.description}
251+ </p>
252+ <table className="data">
253+ <tbody>
254+ <tr>
255+ <th className="left">a</th>
256+ <td>{inst.a}</td>
257+ <th className="left">k</th>
258+ <td>{inst.k}</td>
259+ <th className="left">d</th>
260+ <td>{inst.d}</td>
261+ </tr>
262+ </tbody>
263+ </table>
264+ </div>
265+ <DomainView inst={inst} />
266+ </div>
267+
268+ <h2>Work-precision results</h2>
269+ {committedError && (
270+ <p className="small muted">
271+ Committed results could not be loaded ({committedError}); showing
272+ local runs only.
273+ </p>
274+ )}
275+ <WorkPrecisionChart curves={curves} />
276+ <div className="row" style={{ marginTop: 14, alignItems: "center" }}>
277+ {SOLVERS.map((s) => (
278+ <span key={s.id} style={{ whiteSpace: "nowrap" }}>
279+ <label>
280+ <input
281+ type="checkbox"
282+ checked={!hidden.has(s.id)}
283+ onChange={(e) => {
284+ setHidden((h) => {
285+ const next = new Set(h);
286+ if (e.target.checked) next.delete(s.id);
287+ else next.add(s.id);
288+ return next;
289+ });
290+ }}
291+ />{" "}
292+ <span
293+ className="legend-swatch"
294+ style={{ background: solverColorVar(s.id, allSolverIds) }}
295+ />
296+ {s.name}
297+ </label>{" "}
298+ <button
299+ onClick={() => runSolver(s.id)}
300+ disabled={running !== null}
301+ title={`Run the full ${s.id} sweep on ${instanceId} in this browser`}
302+ >
303+ {running === s.id ? "running…" : "Run in this browser"}
304+ </button>
305+ </span>
306+ ))}
307+ <label>
308+ repeats{" "}
309+ <select
310+ value={repeats}
311+ onChange={(e) => setRepeats(parseInt(e.target.value, 10))}
312+ >
313+ {[1, 3, 5].map((r) => (
314+ <option key={r} value={r}>
315+ {r}
316+ </option>
317+ ))}
318+ </select>
319+ </label>
320+ </div>
321+ {runStatus && <p className="small muted">{runStatus}</p>}
322+ <div className="row" style={{ marginTop: 10, alignItems: "center" }}>
323+ <label className="small">
324+ machine label{" "}
325+ <input
326+ type="text"
327+ placeholder="e.g. office workstation"
328+ value={machineLabel}
329+ onChange={(e) => setMachineLabel(e.target.value)}
330+ />
331+ </label>
332+ {localRuns
333+ .filter((r) => r.done && r.instanceId === instanceId)
334+ .map((r) => (
335+ <button key={r.key} onClick={() => downloadRun(r)}>
336+ Download {r.solverId} result JSON
337+ </button>
338+ ))}
339+ <label className="small">
340+ load result file{" "}
341+ <input
342+ type="file"
343+ accept=".json,application/json"
344+ multiple
345+ onChange={(e) => loadFiles(e.target.files)}
346+ />
347+ </label>
348+ </div>
349+ <PointsTable curves={curves} />
350+
351+ <h2>Solution and error</h2>
352+ <p className="small muted" style={{ maxWidth: 640 }}>
353+ Compute one solve at a chosen resolution and compare the field with
354+ the exact solution. The solution uses a diverging scale about zero;
355+ the error is the absolute pointwise difference on a log scale.
356+ </p>
357+ <SolutionSection inst={inst} />
358+
359+ <h2>Run it outside the browser</h2>
360+ <p style={{ maxWidth: 720 }}>
361+ The same harness runs in node, with the same solvers, the same
362+ protocol, and the same result format (node 20 or newer; no install
363+ step):
364+ </p>
365+ <pre>{`npx ${cliUrl} run --label "my workstation"`}</pre>
366+ <p style={{ maxWidth: 720 }}>
367+ This writes one result JSON per instance and solver. To benchmark
368+ your own solver, point the harness at a MATLAB function file that
369+ implements the problem's solver interface:
370+ </p>
371+ <pre>{`npx ${cliUrl} run --solver-file my_method.m --solver-id my-method`}</pre>
372+ <p style={{ maxWidth: 720 }}>
373+ Result files can be loaded above (load result file) to view them
374+ against the committed curves before submitting anything. To publish,
375+ open a pull request adding the files under <code>results/</code> in
376+ the <a href={RESULTS_REPO_URL}>results repository</a>; provenance
377+ (machine, runtime, numbl version, solver version) travels inside each
378+ file. Solvers in other languages are planned to enter the same way:
379+ run offline, produce the same result format, submit by PR, with the
380+ file marked as not reproducible in the browser.
381+ </p>
382+
383+ <footer>
384+ fastandaccurate · Apache-2.0 ·{" "}
385+ <a href={REPO_URL}>concept-collection/fastandaccurate</a> · numbl{" "}
386+ {__NUMBL_VERSION__}
387+ </footer>
388+ </main>
389+ );
390+}
src/app/colors.tsadded+74−0View file
@@ -0,0 +1,74 @@
1+// Series color assignment and the ramps used by the field views.
2+// Series slots follow the solver identity in fixed order (registry solvers
3+// first, then any loaded solver ids alphabetically); a filter or rerun
4+// never repaints a surviving series.
5+
6+import { SOLVERS } from "../solvers";
7+
8+const SERIES_VARS = [
9+ "--series-1",
10+ "--series-2",
11+ "--series-3",
12+ "--series-4",
13+ "--series-5",
14+];
15+
16+export function solverColorVar(solverId: string, extraIds: string[]): string {
17+ const known = SOLVERS.map((s) => s.id);
18+ const extras = [...new Set(extraIds.filter((id) => !known.includes(id)))].sort();
19+ const order = [...known, ...extras];
20+ const idx = order.indexOf(solverId);
21+ return `var(${SERIES_VARS[Math.max(0, idx) % SERIES_VARS.length]})`;
22+}
23+
24+// Ramps for canvas rendering (canvas cannot read CSS variables per pixel).
25+// Values are the reference palette's sequential blue steps and diverging
26+// blue/red pair, in light- and dark-mode steppings.
27+
28+const SEQ_LIGHT = ["#cde2fb", "#9ec5f4", "#6da7ec", "#3987e5", "#256abf", "#184f95", "#0d366b"];
29+const SEQ_DARK = ["#0d366b", "#184f95", "#256abf", "#3987e5", "#6da7ec", "#9ec5f4", "#cde2fb"];
30+
31+const DIV_LIGHT = { neg: "#2a78d6", mid: "#f0efec", pos: "#e34948" };
32+const DIV_DARK = { neg: "#3987e5", mid: "#383835", pos: "#e66767" };
33+
34+function hexToRgb(hex: string): [number, number, number] {
35+ return [
36+ parseInt(hex.slice(1, 3), 16),
37+ parseInt(hex.slice(3, 5), 16),
38+ parseInt(hex.slice(5, 7), 16),
39+ ];
40+}
41+
42+function lerp(a: number, b: number, t: number): number {
43+ return a + (b - a) * t;
44+}
45+
46+function interpStops(stops: string[], t: number): [number, number, number] {
47+ const x = Math.min(1, Math.max(0, t)) * (stops.length - 1);
48+ const i = Math.min(stops.length - 2, Math.floor(x));
49+ const f = x - i;
50+ const c0 = hexToRgb(stops[i]);
51+ const c1 = hexToRgb(stops[i + 1]);
52+ return [lerp(c0[0], c1[0], f), lerp(c0[1], c1[1], f), lerp(c0[2], c1[2], f)];
53+}
54+
55+export function isDarkMode(): boolean {
56+ return window.matchMedia("(prefers-color-scheme: dark)").matches;
57+}
58+
59+/** Sequential ramp (magnitude), t in [0, 1], light means small. */
60+export function sequentialColor(t: number, dark: boolean): [number, number, number] {
61+ return interpStops(dark ? SEQ_DARK : SEQ_LIGHT, t);
62+}
63+
64+/** Diverging ramp (polarity), t in [-1, 1], gray at 0. */
65+export function divergingColor(t: number, dark: boolean): [number, number, number] {
66+ const d = dark ? DIV_DARK : DIV_LIGHT;
67+ const tt = Math.min(1, Math.max(-1, t));
68+ if (tt < 0) return interpStops([d.neg, d.mid], 1 + tt);
69+ return interpStops([d.mid, d.pos], tt);
70+}
71+
72+export function cssColor(rgb: [number, number, number]): string {
73+ return `rgb(${Math.round(rgb[0])}, ${Math.round(rgb[1])}, ${Math.round(rgb[2])})`;
74+}
src/app/components/DomainView.tsxadded+92−0View file
@@ -0,0 +1,92 @@
1+// Draws the problem geometry: boundary curve, evaluation points, and the
2+// exact solution's sources.
3+
4+import { useEffect, useRef } from "react";
5+import type { Laplace2dInstance } from "../../problems/laplace2d/spec";
6+import {
7+ boundaryPoint,
8+ evalPoints,
9+ sources,
10+} from "../../problems/laplace2d/exact";
11+
12+function token(name: string): string {
13+ return getComputedStyle(document.documentElement).getPropertyValue(name).trim();
14+}
15+
16+export function DomainView({ inst }: { inst: Laplace2dInstance }) {
17+ const canvasRef = useRef<HTMLCanvasElement>(null);
18+
19+ useEffect(() => {
20+ const canvas = canvasRef.current;
21+ if (!canvas) return;
22+ const draw = () => {
23+ const size = 360;
24+ const dpr = window.devicePixelRatio || 1;
25+ canvas.width = size * dpr;
26+ canvas.height = size * dpr;
27+ canvas.style.width = `${size}px`;
28+ canvas.style.height = `${size}px`;
29+ const ctx = canvas.getContext("2d");
30+ if (!ctx) return;
31+ ctx.scale(dpr, dpr);
32+ ctx.clearRect(0, 0, size, size);
33+
34+ const extent = 1 + Math.abs(inst.a) + inst.d + 0.25;
35+ const s = size / (2 * extent);
36+ const X = (x: number) => size / 2 + x * s;
37+ const Y = (y: number) => size / 2 - y * s;
38+
39+ // domain fill + boundary
40+ ctx.beginPath();
41+ for (let i = 0; i <= 512; i++) {
42+ const t = (2 * Math.PI * i) / 512;
43+ const p = boundaryPoint(inst, t);
44+ if (i === 0) ctx.moveTo(X(p.x), Y(p.y));
45+ else ctx.lineTo(X(p.x), Y(p.y));
46+ }
47+ ctx.closePath();
48+ ctx.fillStyle = token("--surface-2");
49+ ctx.fill();
50+ ctx.strokeStyle = token("--text");
51+ ctx.lineWidth = 2;
52+ ctx.stroke();
53+
54+ // evaluation points
55+ ctx.fillStyle = token("--text-2");
56+ for (const p of evalPoints(inst)) {
57+ ctx.beginPath();
58+ ctx.arc(X(p.x), Y(p.y), 2, 0, 2 * Math.PI);
59+ ctx.fill();
60+ }
61+
62+ // sources
63+ ctx.strokeStyle = token("--series-2");
64+ ctx.lineWidth = 2;
65+ for (const src of sources(inst)) {
66+ const cx = X(src.x);
67+ const cy = Y(src.y);
68+ ctx.beginPath();
69+ ctx.moveTo(cx - 5, cy - 5);
70+ ctx.lineTo(cx + 5, cy + 5);
71+ ctx.moveTo(cx - 5, cy + 5);
72+ ctx.lineTo(cx + 5, cy - 5);
73+ ctx.stroke();
74+ }
75+ };
76+ draw();
77+ const mq = window.matchMedia("(prefers-color-scheme: dark)");
78+ mq.addEventListener("change", draw);
79+ return () => mq.removeEventListener("change", draw);
80+ }, [inst]);
81+
82+ return (
83+ <figure style={{ margin: 0 }}>
84+ <canvas ref={canvasRef} />
85+ <figcaption className="field-caption" style={{ maxWidth: 360 }}>
86+ The domain, the 65 evaluation points where solutions are scored
87+ (dots), and the exact solution's sources a distance {inst.d} outside
88+ the boundary (crosses).
89+ </figcaption>
90+ </figure>
91+ );
92+}
src/app/components/FieldView.tsxadded+134−0View file
@@ -0,0 +1,134 @@
1+// Renders a scalar field on the visualization grid as a masked heatmap
2+// with a colorbar. Two modes: "diverging" for signed fields (the solution)
3+// and "logmag" for magnitudes on a log scale (the error).
4+
5+import { useEffect, useRef } from "react";
6+import type { Laplace2dInstance } from "../../problems/laplace2d/spec";
7+import { insideDomain, vizGrid, VIZ_NGRID } from "../../problems/laplace2d/exact";
8+import { cssColor, divergingColor, isDarkMode, sequentialColor } from "../colors";
9+
10+export interface FieldViewProps {
11+ inst: Laplace2dInstance;
12+ /** Values on the viz grid, index p = ix * ngrid + iy. */
13+ values: Float64Array;
14+ mode: "diverging" | "logmag";
15+ title: string;
16+ caption?: string;
17+}
18+
19+export function FieldView({ inst, values, mode, title, caption }: FieldViewProps) {
20+ const canvasRef = useRef<HTMLCanvasElement>(null);
21+ const barRef = useRef<HTMLCanvasElement>(null);
22+ const rangeRef = useRef<HTMLDivElement>(null);
23+
24+ useEffect(() => {
25+ const canvas = canvasRef.current;
26+ const bar = barRef.current;
27+ if (!canvas || !bar) return;
28+
29+ const draw = () => {
30+ const dark = isDarkMode();
31+ const ngrid = VIZ_NGRID;
32+ const { xs } = vizGrid(inst);
33+ // range
34+ let vmin = Infinity;
35+ let vmax = -Infinity;
36+ for (let ix = 0; ix < ngrid; ix++) {
37+ for (let iy = 0; iy < ngrid; iy++) {
38+ if (!insideDomain(inst, xs[ix], xs[iy])) continue;
39+ const v = values[ix * ngrid + iy];
40+ if (!isFinite(v)) continue;
41+ vmin = Math.min(vmin, v);
42+ vmax = Math.max(vmax, v);
43+ }
44+ }
45+ let lo: number;
46+ let hi: number;
47+ let scale: (v: number) => number;
48+ if (mode === "diverging") {
49+ const vabs = Math.max(Math.abs(vmin), Math.abs(vmax), 1e-300);
50+ lo = -vabs;
51+ hi = vabs;
52+ scale = (v) => v / vabs; // [-1, 1]
53+ } else {
54+ hi = Math.max(vmax, 1e-300);
55+ lo = Math.max(vmin, hi * 1e-8, 1e-300);
56+ const llo = Math.log10(lo);
57+ const lhi = Math.log10(hi);
58+ scale = (v) =>
59+ (Math.log10(Math.min(Math.max(v, lo), hi)) - llo) / Math.max(lhi - llo, 1e-12);
60+ }
61+
62+ canvas.width = ngrid;
63+ canvas.height = ngrid;
64+ const ctx = canvas.getContext("2d");
65+ if (!ctx) return;
66+ const img = ctx.createImageData(ngrid, ngrid);
67+ for (let ix = 0; ix < ngrid; ix++) {
68+ for (let iy = 0; iy < ngrid; iy++) {
69+ const px = ix;
70+ const py = ngrid - 1 - iy;
71+ const o = (py * ngrid + px) * 4;
72+ if (!insideDomain(inst, xs[ix], xs[iy])) {
73+ img.data[o + 3] = 0;
74+ continue;
75+ }
76+ const v = values[ix * ngrid + iy];
77+ const rgb =
78+ mode === "diverging"
79+ ? divergingColor(scale(v), dark)
80+ : sequentialColor(scale(v), dark);
81+ img.data[o] = Math.round(rgb[0]);
82+ img.data[o + 1] = Math.round(rgb[1]);
83+ img.data[o + 2] = Math.round(rgb[2]);
84+ img.data[o + 3] = 255;
85+ }
86+ }
87+ ctx.putImageData(img, 0, 0);
88+
89+ // colorbar
90+ const bw = 220;
91+ const bh = 10;
92+ bar.width = bw;
93+ bar.height = bh;
94+ const bctx = bar.getContext("2d");
95+ if (!bctx) return;
96+ for (let i = 0; i < bw; i++) {
97+ const t = i / (bw - 1);
98+ const rgb =
99+ mode === "diverging"
100+ ? divergingColor(2 * t - 1, dark)
101+ : sequentialColor(t, dark);
102+ bctx.fillStyle = cssColor(rgb);
103+ bctx.fillRect(i, 0, 1, bh);
104+ }
105+ if (rangeRef.current) {
106+ const fmt = (v: number) =>
107+ mode === "logmag" ? v.toExponential(1) : v.toPrecision(3);
108+ rangeRef.current.textContent = `${fmt(lo)} … ${fmt(hi)}`;
109+ }
110+ };
111+
112+ draw();
113+ const mq = window.matchMedia("(prefers-color-scheme: dark)");
114+ mq.addEventListener("change", draw);
115+ return () => mq.removeEventListener("change", draw);
116+ }, [inst, values, mode]);
117+
118+ return (
119+ <figure style={{ margin: 0 }}>
120+ <div className="small" style={{ fontWeight: 600, marginBottom: 4 }}>
121+ {title}
122+ </div>
123+ <canvas
124+ ref={canvasRef}
125+ style={{ width: 300, height: 300, imageRendering: "auto" }}
126+ />
127+ <div style={{ display: "flex", alignItems: "center", gap: 8, marginTop: 4 }}>
128+ <canvas ref={barRef} style={{ width: 220, height: 10, borderRadius: 3 }} />
129+ </div>
130+ <div ref={rangeRef} className="small muted" />
131+ {caption && <figcaption className="field-caption" style={{ maxWidth: 300 }}>{caption}</figcaption>}
132+ </figure>
133+ );
134+}
src/app/components/PointsTable.tsxadded+42−0View file
@@ -0,0 +1,42 @@
1+// Table view of the charted data, one collapsible table per curve.
2+
3+import type { ChartCurve } from "./WorkPrecisionChart";
4+
5+export function PointsTable({ curves }: { curves: ChartCurve[] }) {
6+ const nonEmpty = curves.filter((c) => c.points.length > 0);
7+ if (nonEmpty.length === 0) return null;
8+ return (
9+ <div style={{ marginTop: 12 }}>
10+ {nonEmpty.map((c) => (
11+ <details key={c.key} style={{ marginBottom: 6 }}>
12+ <summary className="small" style={{ cursor: "pointer" }}>
13+ <span className="legend-swatch" style={{ background: c.color }} />
14+ {c.label} — data table
15+ </summary>
16+ <table className="data" style={{ marginTop: 6 }}>
17+ <thead>
18+ <tr>
19+ <th>n</th>
20+ <th>rel max error</th>
21+ <th>rel L2 error</th>
22+ <th>solve (s)</th>
23+ <th>cold (s)</th>
24+ </tr>
25+ </thead>
26+ <tbody>
27+ {c.points.map((p) => (
28+ <tr key={p.n}>
29+ <td>{p.n}</td>
30+ <td>{p.relMax.toExponential(2)}</td>
31+ <td>{p.relL2.toExponential(2)}</td>
32+ <td>{p.solveSeconds.toFixed(4)}</td>
33+ <td>{p.coldSeconds.toFixed(4)}</td>
34+ </tr>
35+ ))}
36+ </tbody>
37+ </table>
38+ </details>
39+ ))}
40+ </div>
41+ );
42+}
src/app/components/SolutionSection.tsxadded+132−0View file
@@ -0,0 +1,132 @@
1+// The solution-and-error view: compute one solve in the browser at a
2+// chosen resolution and show the field next to its pointwise error
3+// against the exact solution.
4+
5+import { useMemo, useState } from "react";
6+import type { Laplace2dInstance } from "../../problems/laplace2d/spec";
7+import { exactU, vizGrid, VIZ_NGRID } from "../../problems/laplace2d/exact";
8+import { SOLVERS, getSolver } from "../../solvers";
9+import { solutionInBrowser } from "../workerClient";
10+import type { ResultPoint } from "../../harness/resultSchema";
11+import { FieldView } from "./FieldView";
12+
13+function exactGridValues(inst: Laplace2dInstance): Float64Array {
14+ const { xs } = vizGrid(inst);
15+ const out = new Float64Array(VIZ_NGRID * VIZ_NGRID);
16+ for (let ix = 0; ix < VIZ_NGRID; ix++) {
17+ for (let iy = 0; iy < VIZ_NGRID; iy++) {
18+ out[ix * VIZ_NGRID + iy] = exactU(inst, xs[ix], xs[iy]);
19+ }
20+ }
21+ return out;
22+}
23+
24+interface Computed {
25+ solverId: string;
26+ n: number;
27+ uGrid: Float64Array;
28+ point: ResultPoint;
29+}
30+
31+export function SolutionSection({ inst }: { inst: Laplace2dInstance }) {
32+ const [solverId, setSolverId] = useState(SOLVERS[0].id);
33+ const [n, setN] = useState<number>(
34+ SOLVERS[0].sweepN[Math.floor(SOLVERS[0].sweepN.length * 0.7)]
35+ );
36+ const [busy, setBusy] = useState(false);
37+ const [error, setError] = useState<string | null>(null);
38+ const [computed, setComputed] = useState<Computed | null>(null);
39+
40+ const exact = useMemo(() => exactGridValues(inst), [inst]);
41+
42+ const errField = useMemo(() => {
43+ if (!computed) return null;
44+ const out = new Float64Array(exact.length);
45+ for (let i = 0; i < exact.length; i++) {
46+ out[i] = Math.abs(computed.uGrid[i] - exact[i]);
47+ }
48+ return out;
49+ }, [computed, exact]);
50+
51+ const solver = getSolver(solverId);
52+ const showsComputed = computed !== null && computed.solverId === solverId;
53+
54+ async function compute() {
55+ setBusy(true);
56+ setError(null);
57+ try {
58+ const { point, uGrid } = await solutionInBrowser(inst.id, solverId, n);
59+ setComputed({ solverId, n, uGrid, point });
60+ } catch (err) {
61+ setError(err instanceof Error ? err.message : String(err));
62+ } finally {
63+ setBusy(false);
64+ }
65+ }
66+
67+ return (
68+ <div>
69+ <div className="row" style={{ alignItems: "center", gap: 12, marginBottom: 12 }}>
70+ <label>
71+ solver{" "}
72+ <select
73+ value={solverId}
74+ onChange={(e) => {
75+ const id = e.target.value;
76+ setSolverId(id);
77+ const sw = getSolver(id).sweepN;
78+ setN(sw[Math.floor(sw.length * 0.7)]);
79+ }}
80+ >
81+ {SOLVERS.map((s) => (
82+ <option key={s.id} value={s.id}>
83+ {s.name}
84+ </option>
85+ ))}
86+ </select>
87+ </label>
88+ <label>
89+ n{" "}
90+ <select value={n} onChange={(e) => setN(parseInt(e.target.value, 10))}>
91+ {solver.sweepN.map((v) => (
92+ <option key={v} value={v}>
93+ {v}
94+ </option>
95+ ))}
96+ </select>
97+ </label>
98+ <button className="primary" onClick={compute} disabled={busy}>
99+ {busy ? "computing…" : "Compute in this browser"}
100+ </button>
101+ </div>
102+ {error && <p className="small" style={{ color: "var(--series-2)" }}>{error}</p>}
103+ <div className="row">
104+ <FieldView
105+ inst={inst}
106+ values={exact}
107+ mode="diverging"
108+ title="Exact solution"
109+ caption="The manufactured harmonic function, sampled on the visualization grid."
110+ />
111+ {showsComputed && computed && (
112+ <FieldView
113+ inst={inst}
114+ values={computed.uGrid}
115+ mode="diverging"
116+ title={`${solver.name}, n = ${computed.n}`}
117+ caption={`Computed in this browser: rel max error ${computed.point.relMax.toExponential(2)}, solve ${(computed.point.solveSeconds * 1000).toPrecision(3)} ms.`}
118+ />
119+ )}
120+ {showsComputed && errField && computed && (
121+ <FieldView
122+ inst={inst}
123+ values={errField}
124+ mode="logmag"
125+ title="Pointwise error (log scale)"
126+ caption="Absolute difference from the exact solution; the color scale spans 8 decades below the maximum."
127+ />
128+ )}
129+ </div>
130+ </div>
131+ );
132+}
src/app/components/WorkPrecisionChart.tsxadded+279−0View file
@@ -0,0 +1,279 @@
1+// The work-precision chart: relative max error against median solve time,
2+// both log scale. One curve per (solver, environment); color follows the
3+// solver, line style distinguishes committed / local / loaded origins.
4+
5+import { useMemo, useRef, useState } from "react";
6+import type { ResultPoint } from "../../harness/resultSchema";
7+
8+export interface ChartCurve {
9+ key: string;
10+ solverId: string;
11+ label: string;
12+ color: string;
13+ /** stroke-dasharray, undefined for solid (committed results). */
14+ dash?: string;
15+ /** Open markers (used for runs made in this browser). */
16+ open?: boolean;
17+ points: ResultPoint[];
18+}
19+
20+interface Hover {
21+ px: number;
22+ py: number;
23+ curve: ChartCurve;
24+ point: ResultPoint;
25+}
26+
27+const W = 760;
28+const H = 470;
29+const M = { l: 64, r: 150, t: 14, b: 50 };
30+
31+function decades(min: number, max: number): number[] {
32+ const lo = Math.floor(Math.log10(min));
33+ const hi = Math.ceil(Math.log10(max));
34+ const out: number[] = [];
35+ for (let e = lo; e <= hi; e++) out.push(e);
36+ return out;
37+}
38+
39+function fmtPow(e: number): string {
40+ const sup = String(e)
41+ .split("")
42+ .map(
43+ (c) =>
44+ ({ "-": "⁻", "0": "⁰", "1": "¹", "2": "²", "3": "³", "4": "⁴", "5": "⁵", "6": "⁶", "7": "⁷", "8": "⁸", "9": "⁹" })[c] ?? c
45+ )
46+ .join("");
47+ return `10${sup}`;
48+}
49+
50+export function WorkPrecisionChart({ curves }: { curves: ChartCurve[] }) {
51+ const wrapRef = useRef<HTMLDivElement>(null);
52+ const [hover, setHover] = useState<Hover | null>(null);
53+
54+ const nonEmpty = curves.filter((c) => c.points.length > 0);
55+
56+ const scales = useMemo(() => {
57+ let xMin = Infinity,
58+ xMax = -Infinity,
59+ yMin = Infinity,
60+ yMax = -Infinity;
61+ for (const c of nonEmpty) {
62+ for (const p of c.points) {
63+ const x = Math.max(p.solveSeconds, 1e-8);
64+ const y = Math.max(p.relMax, 1e-17);
65+ xMin = Math.min(xMin, x);
66+ xMax = Math.max(xMax, x);
67+ yMin = Math.min(yMin, y);
68+ yMax = Math.max(yMax, y);
69+ }
70+ }
71+ if (!isFinite(xMin)) {
72+ xMin = 1e-4; xMax = 1; yMin = 1e-12; yMax = 1;
73+ }
74+ const xE = decades(xMin, xMax * 1.0001);
75+ const yE = decades(yMin, yMax * 1.0001);
76+ const xLo = xE[0];
77+ const xHi = xE[xE.length - 1];
78+ const yLo = yE[0];
79+ const yHi = yE[yE.length - 1];
80+ const sx = (v: number) =>
81+ M.l + ((Math.log10(Math.max(v, 1e-17)) - xLo) / Math.max(xHi - xLo, 1)) * (W - M.l - M.r);
82+ const sy = (v: number) =>
83+ H - M.b - ((Math.log10(Math.max(v, 1e-17)) - yLo) / Math.max(yHi - yLo, 1)) * (H - M.t - M.b);
84+ return { sx, sy, xE, yE };
85+ }, [nonEmpty]);
86+
87+ const { sx, sy, xE, yE } = scales;
88+ const xStep = xE.length > 8 ? 2 : 1;
89+ const yStep = yE.length > 8 ? 2 : 1;
90+
91+ function onMove(e: React.MouseEvent<SVGSVGElement>) {
92+ const svg = e.currentTarget;
93+ const rect = svg.getBoundingClientRect();
94+ const px = ((e.clientX - rect.left) / rect.width) * W;
95+ const py = ((e.clientY - rect.top) / rect.height) * H;
96+ let best: Hover | null = null;
97+ let bestD = 26 * 26;
98+ for (const c of nonEmpty) {
99+ for (const p of c.points) {
100+ const dx = sx(p.solveSeconds) - px;
101+ const dy = sy(Math.max(p.relMax, 1e-17)) - py;
102+ const d = dx * dx + dy * dy;
103+ if (d < bestD) {
104+ bestD = d;
105+ best = { px: sx(p.solveSeconds), py: sy(Math.max(p.relMax, 1e-17)), curve: c, point: p };
106+ }
107+ }
108+ }
109+ setHover(best);
110+ }
111+
112+ if (nonEmpty.length === 0) {
113+ return (
114+ <div className="panel muted" style={{ padding: "40px 20px", textAlign: "center" }}>
115+ No results for this selection yet. Run a solver below, or load a
116+ result file produced by the command line.
117+ </div>
118+ );
119+ }
120+
121+ return (
122+ <div ref={wrapRef} style={{ position: "relative", maxWidth: 860 }}>
123+ <svg
124+ viewBox={`0 0 ${W} ${H}`}
125+ style={{ width: "100%", height: "auto", display: "block" }}
126+ onMouseMove={onMove}
127+ onMouseLeave={() => setHover(null)}
128+ role="img"
129+ aria-label="Work-precision chart: relative max error versus median solve time, log-log"
130+ >
131+ {/* grid */}
132+ {xE.map((ex) => (
133+ <line
134+ key={`gx${ex}`}
135+ x1={sx(10 ** ex)}
136+ x2={sx(10 ** ex)}
137+ y1={M.t}
138+ y2={H - M.b}
139+ stroke="var(--grid)"
140+ strokeWidth={1}
141+ />
142+ ))}
143+ {yE.map((ey) => (
144+ <line
145+ key={`gy${ey}`}
146+ x1={M.l}
147+ x2={W - M.r}
148+ y1={sy(10 ** ey)}
149+ y2={sy(10 ** ey)}
150+ stroke="var(--grid)"
151+ strokeWidth={1}
152+ />
153+ ))}
154+ {/* axes */}
155+ <line x1={M.l} x2={W - M.r} y1={H - M.b} y2={H - M.b} stroke="var(--border)" />
156+ <line x1={M.l} x2={M.l} y1={M.t} y2={H - M.b} stroke="var(--border)" />
157+ {xE.map(
158+ (ex, i) =>
159+ i % xStep === 0 && (
160+ <text
161+ key={`tx${ex}`}
162+ x={sx(10 ** ex)}
163+ y={H - M.b + 18}
164+ textAnchor="middle"
165+ fontSize={12}
166+ fill="var(--text-2)"
167+ >
168+ {fmtPow(ex)}
169+ </text>
170+ )
171+ )}
172+ {yE.map(
173+ (ey, i) =>
174+ i % yStep === 0 && (
175+ <text
176+ key={`ty${ey}`}
177+ x={M.l - 8}
178+ y={sy(10 ** ey) + 4}
179+ textAnchor="end"
180+ fontSize={12}
181+ fill="var(--text-2)"
182+ >
183+ {fmtPow(ey)}
184+ </text>
185+ )
186+ )}
187+ <text
188+ x={(M.l + W - M.r) / 2}
189+ y={H - 10}
190+ textAnchor="middle"
191+ fontSize={13}
192+ fill="var(--text-2)"
193+ >
194+ median solve time (seconds)
195+ </text>
196+ <text
197+ x={16}
198+ y={(M.t + H - M.b) / 2}
199+ textAnchor="middle"
200+ fontSize={13}
201+ fill="var(--text-2)"
202+ transform={`rotate(-90 16 ${(M.t + H - M.b) / 2})`}
203+ >
204+ relative max error
205+ </text>
206+ {/* curves */}
207+ {nonEmpty.map((c) => {
208+ const pts = [...c.points].sort((a, b) => a.solveSeconds - b.solveSeconds);
209+ const path = pts
210+ .map(
211+ (p, i) =>
212+ `${i === 0 ? "M" : "L"}${sx(p.solveSeconds).toFixed(1)},${sy(Math.max(p.relMax, 1e-17)).toFixed(1)}`
213+ )
214+ .join("");
215+ const last = pts[pts.length - 1];
216+ return (
217+ <g key={c.key}>
218+ <path
219+ d={path}
220+ fill="none"
221+ stroke={c.color}
222+ strokeWidth={2}
223+ strokeDasharray={c.dash}
224+ />
225+ {pts.map((p, i) => (
226+ <circle
227+ key={i}
228+ cx={sx(p.solveSeconds)}
229+ cy={sy(Math.max(p.relMax, 1e-17))}
230+ r={4}
231+ fill={c.open ? "var(--surface)" : c.color}
232+ stroke={c.color}
233+ strokeWidth={c.open ? 2 : 0}
234+ />
235+ ))}
236+ {nonEmpty.length <= 4 && (
237+ <text
238+ x={sx(last.solveSeconds) + 9}
239+ y={sy(Math.max(last.relMax, 1e-17)) + 4}
240+ fontSize={12}
241+ fill="var(--text-2)"
242+ >
243+ {c.label}
244+ </text>
245+ )}
246+ </g>
247+ );
248+ })}
249+ {hover && (
250+ <circle
251+ cx={hover.px}
252+ cy={hover.py}
253+ r={7}
254+ fill="none"
255+ stroke={hover.curve.color}
256+ strokeWidth={2}
257+ />
258+ )}
259+ </svg>
260+ {hover && wrapRef.current && (
261+ <div
262+ className="chart-tooltip"
263+ style={{
264+ left: `${(hover.px / W) * wrapRef.current.clientWidth + 12}px`,
265+ top: `${(hover.py / H) * (wrapRef.current.clientWidth * (H / W)) - 10}px`,
266+ }}
267+ >
268+ <div>
269+ <strong>{hover.curve.label}</strong>
270+ </div>
271+ <div>n = {hover.point.n}</div>
272+ <div>rel max error = {hover.point.relMax.toExponential(2)}</div>
273+ <div>rel L2 error = {hover.point.relL2.toExponential(2)}</div>
274+ <div>solve = {(hover.point.solveSeconds * 1000).toPrecision(3)} ms</div>
275+ </div>
276+ )}
277+ </div>
278+ );
279+}
src/app/main.tsxadded+10−0View file
@@ -0,0 +1,10 @@
1+import { StrictMode } from "react";
2+import { createRoot } from "react-dom/client";
3+import { App } from "./App";
4+import "./styles.css";
5+
6+createRoot(document.getElementById("root")!).render(
7+ <StrictMode>
8+ <App />
9+ </StrictMode>
10+);
src/app/matlabSources.tsadded+25−0View file
@@ -0,0 +1,25 @@
1+// The MATLAB sources, inlined into the bundle by vite. Used by the worker
2+// (and only there; the node CLI reads the same files from disk).
3+
4+const files = import.meta.glob("../{problems,solvers}/**/*.m", {
5+ query: "?raw",
6+ import: "default",
7+ eager: true,
8+}) as Record<string, string>;
9+
10+function get(path: string): string {
11+ const src = files[path];
12+ if (!src) throw new Error(`missing MATLAB source: ${path}`);
13+ return src;
14+}
15+
16+export function matlabBase() {
17+ return {
18+ buildProblem: get("../problems/laplace2d/matlab/build_problem.m"),
19+ bdata: get("../problems/laplace2d/matlab/laplace2d_bdata.m"),
20+ };
21+}
22+
23+export function solverSource(solverId: string): string {
24+ return get(`../solvers/${solverId}/solver.m`);
25+}
src/app/results.tsadded+50−0View file
@@ -0,0 +1,50 @@
1+// Committed results live in the fastandaccurate-results repository and are
2+// fetched statically: index.json lists the result files, each of which is
3+// one work-precision sweep in the format of src/harness/resultSchema.ts.
4+
5+import {
6+ RESULT_FORMAT,
7+ type ResultFile,
8+} from "../harness/resultSchema";
9+
10+export const RESULTS_REPO_URL =
11+ "https://github.com/concept-collection/fastandaccurate-results";
12+const RESULTS_RAW =
13+ "https://raw.githubusercontent.com/concept-collection/fastandaccurate-results/main";
14+
15+export function isResultFile(x: unknown): x is ResultFile {
16+ const r = x as ResultFile;
17+ return (
18+ !!r &&
19+ r.format === RESULT_FORMAT &&
20+ typeof r.problem === "string" &&
21+ typeof r.instance === "string" &&
22+ !!r.solver &&
23+ typeof r.solver.id === "string" &&
24+ Array.isArray(r.points)
25+ );
26+}
27+
28+export async function fetchCommittedResults(): Promise<ResultFile[]> {
29+ const idxResp = await fetch(`${RESULTS_RAW}/index.json`, { cache: "no-cache" });
30+ if (!idxResp.ok) throw new Error(`index.json: HTTP ${idxResp.status}`);
31+ const idx = (await idxResp.json()) as { files?: string[] };
32+ const files = idx.files ?? [];
33+ const results = await Promise.all(
34+ files.map(async (f) => {
35+ const resp = await fetch(`${RESULTS_RAW}/${f}`, { cache: "no-cache" });
36+ if (!resp.ok) return null;
37+ const data: unknown = await resp.json();
38+ return isResultFile(data) ? data : null;
39+ })
40+ );
41+ return results.filter((r): r is ResultFile => r !== null);
42+}
43+
44+/** Short human label for the environment a result was measured in. */
45+export function environmentLabel(r: ResultFile): string {
46+ const env = r.environment;
47+ if (env.machineLabel) return `${env.machineLabel} (${env.kind})`;
48+ if (env.kind === "node") return `${env.cpu ?? "unknown cpu"} (node)`;
49+ return "browser";
50+}
src/app/styles.cssadded+230−0View file
@@ -0,0 +1,230 @@
1+/* Design tokens. The series and ramp values are the validated reference
2+ * palette from the dataviz method; series slots are assigned to solvers
3+ * in fixed registry order and never cycled. */
4+
5+:root {
6+ color-scheme: light;
7+ --surface: #fcfcfb;
8+ --surface-2: #f4f3f0;
9+ --text: #0b0b0b;
10+ --text-2: #52514e;
11+ --text-3: #767570;
12+ --border: #dcdbd6;
13+ --grid: #eceae6;
14+ --link: #1c5cab;
15+ --series-1: #2a78d6;
16+ --series-2: #eb6834;
17+ --series-3: #1baf7a;
18+ --series-4: #eda100;
19+ --series-5: #e87ba4;
20+ --accent: #2a78d6;
21+}
22+
23+@media (prefers-color-scheme: dark) {
24+ :root {
25+ color-scheme: dark;
26+ --surface: #1a1a19;
27+ --surface-2: #232322;
28+ --text: #ffffff;
29+ --text-2: #c3c2b7;
30+ --text-3: #8f8e86;
31+ --border: #3a3a37;
32+ --grid: #2c2c2a;
33+ --link: #86b6ef;
34+ --series-1: #3987e5;
35+ --series-2: #d95926;
36+ --series-3: #199e70;
37+ --series-4: #c98500;
38+ --series-5: #d55181;
39+ --accent: #3987e5;
40+ }
41+}
42+
43+* {
44+ box-sizing: border-box;
45+}
46+
47+body {
48+ margin: 0;
49+ background: var(--surface);
50+ color: var(--text);
51+ font-family: system-ui, -apple-system, "Segoe UI", Roboto, sans-serif;
52+ font-size: 15px;
53+ line-height: 1.55;
54+}
55+
56+main {
57+ max-width: 1100px;
58+ margin: 0 auto;
59+ padding: 0 20px 80px;
60+}
61+
62+a {
63+ color: var(--link);
64+}
65+
66+h1 {
67+ font-size: 1.7rem;
68+ font-weight: 650;
69+ margin: 1.6rem 0 0.2rem;
70+}
71+
72+h2 {
73+ font-size: 1.25rem;
74+ font-weight: 650;
75+ margin: 2.6rem 0 0.6rem;
76+ padding-top: 1.2rem;
77+ border-top: 1px solid var(--border);
78+}
79+
80+h3 {
81+ font-size: 1.02rem;
82+ font-weight: 650;
83+ margin: 1.4rem 0 0.4rem;
84+}
85+
86+p {
87+ margin: 0.5rem 0;
88+ max-width: 72ch;
89+}
90+
91+.subtitle {
92+ color: var(--text-2);
93+ margin-top: 0;
94+}
95+
96+code,
97+pre {
98+ font-family: ui-monospace, "SF Mono", Menlo, Consolas, monospace;
99+ font-size: 0.88em;
100+}
101+
102+pre {
103+ background: var(--surface-2);
104+ border: 1px solid var(--border);
105+ border-radius: 6px;
106+ padding: 10px 14px;
107+ overflow-x: auto;
108+}
109+
110+.muted {
111+ color: var(--text-2);
112+}
113+
114+.small {
115+ font-size: 0.85rem;
116+}
117+
118+.row {
119+ display: flex;
120+ flex-wrap: wrap;
121+ gap: 20px;
122+ align-items: flex-start;
123+}
124+
125+.panel {
126+ border: 1px solid var(--border);
127+ border-radius: 8px;
128+ padding: 14px 16px;
129+ background: var(--surface);
130+}
131+
132+button {
133+ font: inherit;
134+ font-size: 0.9rem;
135+ padding: 5px 14px;
136+ border: 1px solid var(--border);
137+ border-radius: 6px;
138+ background: var(--surface-2);
139+ color: var(--text);
140+ cursor: pointer;
141+}
142+
143+button:hover:not(:disabled) {
144+ border-color: var(--text-3);
145+}
146+
147+button:disabled {
148+ opacity: 0.5;
149+ cursor: default;
150+}
151+
152+button.primary {
153+ background: var(--accent);
154+ border-color: var(--accent);
155+ color: #ffffff;
156+}
157+
158+select,
159+input[type="text"] {
160+ font: inherit;
161+ font-size: 0.9rem;
162+ padding: 4px 8px;
163+ border: 1px solid var(--border);
164+ border-radius: 6px;
165+ background: var(--surface);
166+ color: var(--text);
167+}
168+
169+label {
170+ font-size: 0.9rem;
171+}
172+
173+table.data {
174+ border-collapse: collapse;
175+ font-size: 0.85rem;
176+ font-variant-numeric: tabular-nums;
177+}
178+
179+table.data th,
180+table.data td {
181+ border: 1px solid var(--border);
182+ padding: 3px 10px;
183+ text-align: right;
184+}
185+
186+table.data th {
187+ background: var(--surface-2);
188+ font-weight: 600;
189+}
190+
191+table.data td.left,
192+table.data th.left {
193+ text-align: left;
194+}
195+
196+.legend-swatch {
197+ display: inline-block;
198+ width: 14px;
199+ height: 3px;
200+ border-radius: 2px;
201+ vertical-align: middle;
202+ margin-right: 6px;
203+}
204+
205+.chart-tooltip {
206+ position: absolute;
207+ pointer-events: none;
208+ background: var(--surface);
209+ border: 1px solid var(--border);
210+ border-radius: 6px;
211+ padding: 6px 10px;
212+ font-size: 0.8rem;
213+ box-shadow: 0 2px 8px rgba(0, 0, 0, 0.12);
214+ z-index: 10;
215+ white-space: nowrap;
216+}
217+
218+.field-caption {
219+ font-size: 0.85rem;
220+ color: var(--text-2);
221+ margin-top: 4px;
222+}
223+
224+footer {
225+ margin-top: 4rem;
226+ padding-top: 1rem;
227+ border-top: 1px solid var(--border);
228+ color: var(--text-3);
229+ font-size: 0.85rem;
230+}
src/app/vite-env.d.tsadded+4−0View file
@@ -0,0 +1,4 @@
1+/// <reference types="vite/client" />
2+
3+declare const __NUMBL_VERSION__: string;
4+declare const __BUILD_ID__: string;
src/app/worker.tsadded+92−0View file
@@ -0,0 +1,92 @@
1+/// <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.
4+
5+import { getInstance } from "../problems/laplace2d/spec";
6+import { getSolver } from "../solvers";
7+import { runSweep } from "../harness/sweep";
8+import { runPoint } from "../harness/runner";
9+import { toResultPoint, type ResultPoint } from "../harness/resultSchema";
10+import { matlabBase, solverSource } from "./matlabSources";
11+
12+export interface SweepRequest {
13+ type: "sweep";
14+ id: number;
15+ instanceId: string;
16+ solverId: string;
17+ repeats: number;
18+}
19+
20+export interface SolutionRequest {
21+ type: "solution";
22+ id: number;
23+ instanceId: string;
24+ solverId: string;
25+ n: number;
26+}
27+
28+export type WorkerRequest = SweepRequest | SolutionRequest;
29+
30+export type WorkerResponse =
31+ | { type: "point"; id: number; point: ResultPoint; index: number; total: number }
32+ | { type: "sweepDone"; id: number; points: ResultPoint[] }
33+ | {
34+ type: "solutionDone";
35+ id: number;
36+ point: ResultPoint;
37+ uGrid: Float64Array;
38+ }
39+ | { type: "error"; id: number; message: string };
40+
41+self.onmessage = (e: MessageEvent<WorkerRequest>) => {
42+ const msg = e.data;
43+ try {
44+ if (msg.type === "sweep") {
45+ const points = runSweep({
46+ instance: getInstance(msg.instanceId),
47+ solver: getSolver(msg.solverId),
48+ sources: { ...matlabBase(), solver: solverSource(msg.solverId) },
49+ repeats: msg.repeats,
50+ onPoint: (p, index, total) => {
51+ const resp: WorkerResponse = {
52+ type: "point",
53+ id: msg.id,
54+ point: toResultPoint(p),
55+ index,
56+ total,
57+ };
58+ postMessage(resp);
59+ },
60+ });
61+ const resp: WorkerResponse = {
62+ type: "sweepDone",
63+ id: msg.id,
64+ points: points.map(toResultPoint),
65+ };
66+ postMessage(resp);
67+ } else if (msg.type === "solution") {
68+ const p = runPoint({
69+ instance: getInstance(msg.instanceId),
70+ n: msg.n,
71+ repeats: 1,
72+ wantGrid: true,
73+ sources: { ...matlabBase(), solver: solverSource(msg.solverId) },
74+ });
75+ if (!p.uGrid) throw new Error("solver returned no grid values");
76+ const resp: WorkerResponse = {
77+ type: "solutionDone",
78+ id: msg.id,
79+ point: toResultPoint(p),
80+ uGrid: p.uGrid,
81+ };
82+ postMessage(resp, { transfer: [p.uGrid.buffer] });
83+ }
84+ } catch (err) {
85+ const resp: WorkerResponse = {
86+ type: "error",
87+ id: msg.id,
88+ message: err instanceof Error ? err.message : String(err),
89+ };
90+ postMessage(resp);
91+ }
92+};
src/app/workerClient.tsadded+73−0View file
@@ -0,0 +1,73 @@
1+// Host-side wrapper for the compute worker: one shared worker, one
2+// promise per request. Requests are serviced one at a time in order.
3+
4+import type { WorkerRequest, WorkerResponse } from "./worker";
5+import type { ResultPoint } from "../harness/resultSchema";
6+
7+type Pending = {
8+ onPoint?: (point: ResultPoint, index: number, total: number) => void;
9+ resolve: (value: never) => void;
10+ reject: (err: Error) => void;
11+};
12+
13+let worker: Worker | null = null;
14+let nextId = 1;
15+const pending = new Map<number, Pending>();
16+
17+function getWorker(): Worker {
18+ if (!worker) {
19+ worker = new Worker(new URL("./worker.ts", import.meta.url), {
20+ type: "module",
21+ });
22+ worker.onmessage = (e: MessageEvent<WorkerResponse>) => {
23+ const msg = e.data;
24+ const p = pending.get(msg.id);
25+ if (!p) return;
26+ if (msg.type === "point") {
27+ p.onPoint?.(msg.point, msg.index, msg.total);
28+ } else if (msg.type === "sweepDone") {
29+ pending.delete(msg.id);
30+ (p.resolve as (v: ResultPoint[]) => void)(msg.points);
31+ } else if (msg.type === "solutionDone") {
32+ pending.delete(msg.id);
33+ (p.resolve as (v: { point: ResultPoint; uGrid: Float64Array }) => void)({
34+ point: msg.point,
35+ uGrid: msg.uGrid,
36+ });
37+ } else if (msg.type === "error") {
38+ pending.delete(msg.id);
39+ p.reject(new Error(msg.message));
40+ }
41+ };
42+ }
43+ return worker;
44+}
45+
46+function post(req: WorkerRequest) {
47+ getWorker().postMessage(req);
48+}
49+
50+export function sweepInBrowser(
51+ instanceId: string,
52+ solverId: string,
53+ repeats: number,
54+ onPoint: (point: ResultPoint, index: number, total: number) => void
55+): Promise<ResultPoint[]> {
56+ const id = nextId++;
57+ return new Promise((resolve, reject) => {
58+ pending.set(id, { onPoint, resolve: resolve as never, reject });
59+ post({ type: "sweep", id, instanceId, solverId, repeats });
60+ });
61+}
62+
63+export function solutionInBrowser(
64+ instanceId: string,
65+ solverId: string,
66+ n: number
67+): Promise<{ point: ResultPoint; uGrid: Float64Array }> {
68+ const id = nextId++;
69+ return new Promise((resolve, reject) => {
70+ pending.set(id, { resolve: resolve as never, reject });
71+ post({ type: "solution", id, instanceId, solverId, n });
72+ });
73+}
src/cli/main.tsadded+221−0View file
@@ -0,0 +1,221 @@
1+// The fastandaccurate command line: run work-precision sweeps outside the
2+// browser, through the same harness the site uses, and write result JSON
3+// files ready to submit to the fastandaccurate-results repository by PR.
4+//
5+// fastandaccurate list
6+// fastandaccurate run [--instance <id>] [--solver <id>]
7+// [--solver-file f.m --solver-id name]
8+// [--repeats N] [--max-n N] [--label "text"]
9+// [--out dir]
10+//
11+// In development: npx tsx src/cli/main.ts run ...
12+
13+import { readFileSync, writeFileSync, mkdirSync, existsSync } from "fs";
14+import { dirname, join, resolve } from "path";
15+import { fileURLToPath } from "url";
16+import os from "os";
17+import { INSTANCES, getInstance } from "../problems/laplace2d/spec";
18+import { SOLVERS, getSolver, type SolverManifest } from "../solvers";
19+import { runSweep } from "../harness/sweep";
20+import type { MatlabSources } from "../harness/runner";
21+import {
22+ buildResultFile,
23+ toResultPoint,
24+ type ResultEnvironment,
25+} from "../harness/resultSchema";
26+
27+const moduleDir = dirname(fileURLToPath(import.meta.url));
28+
29+/** Locate the directory holding the .m sources: the repo's src/ in
30+ * development, the bundle's own src/ in the packed CLI. */
31+function findSrcRoot(): string {
32+ const candidates = [join(moduleDir, "..", ".."), join(moduleDir)];
33+ for (const c of candidates) {
34+ if (existsSync(join(c, "src", "problems", "laplace2d", "matlab", "build_problem.m"))) {
35+ return join(c, "src");
36+ }
37+ }
38+ throw new Error("cannot locate MATLAB sources next to the CLI");
39+}
40+
41+const srcRoot = findSrcRoot();
42+const readSrc = (p: string) => readFileSync(join(srcRoot, p), "utf-8");
43+
44+function numblVersion(): string {
45+ // The packed CLI has numbl bundled in; the version is stamped at build
46+ // time. In development (tsx), read it from node_modules instead.
47+ if (typeof __NUMBL_VERSION__ !== "undefined") return __NUMBL_VERSION__;
48+ const c = join(srcRoot, "..", "node_modules", "numbl", "package.json");
49+ if (existsSync(c)) {
50+ return (JSON.parse(readFileSync(c, "utf-8")) as { version: string }).version;
51+ }
52+ return "unknown";
53+}
54+
55+function environment(machineLabel: string | undefined, builtin: boolean): ResultEnvironment {
56+ return {
57+ kind: "node",
58+ runtime: `node ${process.version}`,
59+ numblVersion: numblVersion(),
60+ os: `${os.platform()} ${os.release()}`,
61+ cpu: os.cpus()[0]?.model?.trim() ?? "unknown",
62+ machineLabel,
63+ browserReproducible: builtin,
64+ };
65+}
66+
67+interface Args {
68+ command: string;
69+ flags: Record<string, string>;
70+}
71+
72+function parseArgs(argv: string[]): Args {
73+ const [command = "help", ...rest] = argv;
74+ const flags: Record<string, string> = {};
75+ for (let i = 0; i < rest.length; i++) {
76+ const a = rest[i];
77+ if (!a.startsWith("--")) throw new Error(`unexpected argument: ${a}`);
78+ const key = a.slice(2);
79+ const val = rest[i + 1];
80+ if (val === undefined || val.startsWith("--")) {
81+ flags[key] = "true";
82+ } else {
83+ flags[key] = val;
84+ i++;
85+ }
86+ }
87+ return { command, flags };
88+}
89+
90+function listCommand() {
91+ console.log("Problem: laplace-dirichlet-2d (v1)\n");
92+ console.log("Instances:");
93+ for (const inst of INSTANCES) {
94+ console.log(` ${inst.id.padEnd(14)} ${inst.label}`);
95+ }
96+ console.log("\nSolvers:");
97+ for (const s of SOLVERS) {
98+ console.log(` ${s.id.padEnd(14)} ${s.name} (v${s.version}, ${s.backend})`);
99+ }
100+}
101+
102+async function runCommand(flags: Record<string, string>) {
103+ const repeats = flags.repeats ? parseInt(flags.repeats, 10) : 3;
104+ const maxN = flags["max-n"] ? parseInt(flags["max-n"], 10) : undefined;
105+ const outDir = resolve(flags.out ?? "fastandaccurate-results-out");
106+ const instances = flags.instance
107+ ? [getInstance(flags.instance)]
108+ : INSTANCES;
109+
110+ const base = {
111+ buildProblem: readSrc("problems/laplace2d/matlab/build_problem.m"),
112+ bdata: readSrc("problems/laplace2d/matlab/laplace2d_bdata.m"),
113+ };
114+
115+ let solverList: { manifest: SolverManifest; sources: MatlabSources; source: string }[];
116+ if (flags["solver-file"]) {
117+ const file = resolve(flags["solver-file"]);
118+ const id = flags["solver-id"];
119+ if (!id) throw new Error("--solver-file requires --solver-id");
120+ const manifest: SolverManifest = {
121+ id,
122+ name: id,
123+ description: `custom solver from ${file}`,
124+ version: flags["solver-version"] ?? "0.0.0",
125+ backend: "cpu",
126+ sweepN: getSolver("nystrom-dlp").sweepN,
127+ };
128+ solverList = [
129+ {
130+ manifest,
131+ sources: { ...base, solver: readFileSync(file, "utf-8") },
132+ source: file,
133+ },
134+ ];
135+ } else {
136+ const wanted = flags.solver ? [getSolver(flags.solver)] : SOLVERS;
137+ solverList = wanted.map((manifest) => ({
138+ manifest,
139+ sources: { ...base, solver: readSrc(`solvers/${manifest.id}/solver.m`) },
140+ source: "builtin",
141+ }));
142+ }
143+
144+ mkdirSync(outDir, { recursive: true });
145+ const env = environment(flags.label, !flags["solver-file"]);
146+
147+ for (const inst of instances) {
148+ for (const { manifest, sources, source } of solverList) {
149+ console.log(`\n${inst.id} / ${manifest.id}`);
150+ console.log(" n relMax relL2 solve(s)");
151+ const points = runSweep({
152+ instance: inst,
153+ solver: manifest,
154+ sources,
155+ repeats,
156+ maxN,
157+ onPoint: (p) => {
158+ console.log(
159+ ` ${String(p.n).padStart(4)} ${p.relMax.toExponential(2)} ` +
160+ `${p.relL2.toExponential(2)} ${p.solveSeconds.toFixed(4)}`
161+ );
162+ },
163+ });
164+ const result = await buildResultFile({
165+ instance: inst,
166+ solver: {
167+ id: manifest.id,
168+ version: manifest.version,
169+ backend: manifest.backend,
170+ source,
171+ },
172+ environment: env,
173+ repeats,
174+ points: points.map(toResultPoint),
175+ });
176+ const name = `laplace-dirichlet-2d.${inst.id}.${manifest.id}.json`;
177+ const path = join(outDir, name);
178+ writeFileSync(path, JSON.stringify(result, null, 2) + "\n");
179+ console.log(` wrote ${path}`);
180+ }
181+ }
182+ console.log(
183+ "\nTo publish: open a pull request adding these files under results/ in " +
184+ "https://github.com/concept-collection/fastandaccurate-results"
185+ );
186+}
187+
188+async function main() {
189+ const args = parseArgs(process.argv.slice(2));
190+ if (args.command === "list") {
191+ listCommand();
192+ } else if (args.command === "run") {
193+ await runCommand(args.flags);
194+ } else {
195+ console.log(
196+ [
197+ "fastandaccurate - PDE solver benchmarks (https://concept-collection.github.io/fastandaccurate/)",
198+ "",
199+ "Commands:",
200+ " list List problems, instances, and solvers",
201+ " run Run work-precision sweeps and write result JSON files",
202+ "",
203+ "Run flags:",
204+ " --instance <id> One instance (default: all)",
205+ " --solver <id> One built-in solver (default: all)",
206+ " --solver-file <f.m> A custom solver file (requires --solver-id)",
207+ " --solver-id <name> Identifier for the custom solver",
208+ " --solver-version <v> Version string for the custom solver",
209+ " --repeats <N> Timed repeats per point (default 3)",
210+ " --max-n <N> Restrict the sweep to n <= N",
211+ " --label <text> Free-text machine label recorded in results",
212+ " --out <dir> Output directory (default fastandaccurate-results-out)",
213+ ].join("\n")
214+ );
215+ }
216+}
217+
218+main().catch((err) => {
219+ console.error(err instanceof Error ? err.message : err);
220+ process.exit(1);
221+});
src/harness/numblRun.tsadded+65−0View file
@@ -0,0 +1,65 @@
1+// Thin wrapper around numbl's synchronous executeCode for harness runs.
2+// Works identically in a browser worker and in node: no file I/O adapters
3+// are attached, so the MATLAB side must communicate through workspace
4+// variables, which we read back from result.variableValues.
5+
6+import { executeCode } from "numbl";
7+
8+const PROJ = "/fastandaccurate";
9+
10+export interface NumblRunResult {
11+ /** Console output of the run. */
12+ output: string;
13+ /** Named numeric results pulled from the final workspace. */
14+ vars: Record<string, Float64Array>;
15+}
16+
17+/**
18+ * Run mainSource as the main script with the given auxiliary .m files on
19+ * the search path, and extract the requested workspace variables, which
20+ * must be real numeric arrays (or scalars, returned as length-1 arrays).
21+ * Throws on MATLAB errors and on missing/non-numeric variables.
22+ */
23+export function runNumblScript(
24+ mainSource: string,
25+ files: Record<string, string>,
26+ wantVars: string[]
27+): NumblRunResult {
28+ const workspaceFiles = Object.entries(files).map(([name, source]) => ({
29+ name: `${PROJ}/${name}`,
30+ source,
31+ }));
32+ const outputs: string[] = [];
33+ const result = executeCode(
34+ mainSource,
35+ {
36+ onOutput: (text) => outputs.push(text),
37+ displayResults: false,
38+ optimization: "1",
39+ implicitCwdPath: null,
40+ },
41+ workspaceFiles,
42+ `${PROJ}/main.m`,
43+ [PROJ]
44+ );
45+ const vars: Record<string, Float64Array> = {};
46+ for (const name of wantVars) {
47+ const v = result.variableValues[name];
48+ if (typeof v === "number") {
49+ vars[name] = new Float64Array([v]);
50+ } else if (
51+ v &&
52+ typeof v === "object" &&
53+ (v as { kind?: string }).kind === "tensor"
54+ ) {
55+ const tensor = v as { data: Float64Array; imag?: Float64Array };
56+ if (tensor.imag) {
57+ throw new Error(`variable ${name} is complex; expected real`);
58+ }
59+ vars[name] = tensor.data;
60+ } else {
61+ throw new Error(`variable ${name} missing or not numeric after run`);
62+ }
63+ }
64+ return { output: outputs.join(""), vars };
65+}
src/harness/resultSchema.tsadded+112−0View file
@@ -0,0 +1,112 @@
1+// The result-file format committed to the fastandaccurate-results repo.
2+// One file holds one work-precision sweep: one problem instance, one
3+// solver, one environment.
4+
5+import {
6+ canonicalSpec,
7+ canonicalSpecJson,
8+ PROBLEM_ID,
9+ PROBLEM_VERSION,
10+ type Laplace2dInstance,
11+} from "../problems/laplace2d/spec";
12+import type { RunPoint } from "./runner";
13+
14+export const RESULT_FORMAT = "fastandaccurate-result";
15+export const RESULT_FORMAT_VERSION = 1;
16+
17+export interface ResultEnvironment {
18+ kind: "browser" | "node";
19+ /** User agent string (browser) or node version (node). */
20+ runtime: string;
21+ numblVersion: string;
22+ os?: string;
23+ cpu?: string;
24+ /** Free-text label a human recognizes ("office workstation"). */
25+ machineLabel?: string;
26+ /** Whether a visitor can rerun this result in the browser. */
27+ browserReproducible: boolean;
28+}
29+
30+export interface ResultPoint {
31+ n: number;
32+ solveSeconds: number;
33+ solveSecondsAll: number[];
34+ coldSeconds: number;
35+ relMax: number;
36+ relL2: number;
37+}
38+
39+export interface ResultFile {
40+ format: typeof RESULT_FORMAT;
41+ formatVersion: typeof RESULT_FORMAT_VERSION;
42+ problem: string;
43+ problemVersion: number;
44+ instance: string;
45+ spec: Record<string, unknown>;
46+ specHash: string;
47+ solver: {
48+ id: string;
49+ version: string;
50+ backend: "cpu" | "gpu";
51+ /** "builtin" for solvers in the fastandaccurate repo; otherwise a
52+ * URL or free-text pointer to the solver's source. */
53+ source: string;
54+ };
55+ environment: ResultEnvironment;
56+ protocol: {
57+ warmupRuns: number;
58+ timedRuns: number;
59+ timer: string;
60+ };
61+ createdUtc: string;
62+ points: ResultPoint[];
63+}
64+
65+/** SHA-256 hex digest, using WebCrypto (browser, worker, node >= 20). */
66+export async function sha256Hex(text: string): Promise<string> {
67+ const digest = await crypto.subtle.digest(
68+ "SHA-256",
69+ new TextEncoder().encode(text)
70+ );
71+ return Array.from(new Uint8Array(digest))
72+ .map((b) => b.toString(16).padStart(2, "0"))
73+ .join("");
74+}
75+
76+export function toResultPoint(p: RunPoint): ResultPoint {
77+ return {
78+ n: p.n,
79+ solveSeconds: p.solveSeconds,
80+ solveSecondsAll: p.solveSecondsAll,
81+ coldSeconds: p.coldSeconds,
82+ relMax: p.relMax,
83+ relL2: p.relL2,
84+ };
85+}
86+
87+export async function buildResultFile(opts: {
88+ instance: Laplace2dInstance;
89+ solver: { id: string; version: string; backend: "cpu" | "gpu"; source: string };
90+ environment: ResultEnvironment;
91+ repeats: number;
92+ points: ResultPoint[];
93+}): Promise<ResultFile> {
94+ return {
95+ format: RESULT_FORMAT,
96+ formatVersion: RESULT_FORMAT_VERSION,
97+ problem: PROBLEM_ID,
98+ problemVersion: PROBLEM_VERSION,
99+ instance: opts.instance.id,
100+ spec: canonicalSpec(opts.instance) as unknown as Record<string, unknown>,
101+ specHash: await sha256Hex(canonicalSpecJson(opts.instance)),
102+ solver: opts.solver,
103+ environment: opts.environment,
104+ protocol: {
105+ warmupRuns: 1,
106+ timedRuns: opts.repeats,
107+ timer: "numbl tic/toc",
108+ },
109+ createdUtc: new Date().toISOString(),
110+ points: opts.points,
111+ };
112+}
src/harness/runner.tsadded+89−0View file
@@ -0,0 +1,89 @@
1+// Runs one solver at one resolution on one instance and measures it.
2+// The timing protocol (see docs/problems/laplace-dirichlet-2d.md): one
3+// untimed warmup run absorbs JIT compilation, then `repeats` timed runs
4+// whose median is the reported solve time. The warmup itself is timed
5+// and reported as the cold time. All timing is MATLAB tic/toc inside the
6+// numbl session, so browser and node measure the same thing.
7+
8+import { runNumblScript } from "./numblRun";
9+import type { Laplace2dInstance } from "../problems/laplace2d/spec";
10+import { evalErrors } from "../problems/laplace2d/exact";
11+
12+export interface MatlabSources {
13+ /** build_problem.m source */
14+ buildProblem: string;
15+ /** laplace2d_bdata.m source */
16+ bdata: string;
17+ /** solver.m source of the solver under test */
18+ solver: string;
19+}
20+
21+export interface RunPointRequest {
22+ instance: Laplace2dInstance;
23+ n: number;
24+ /** Timed repeats after the warmup (default 3). */
25+ repeats?: number;
26+ /** Also evaluate the solution on the visualization grid. */
27+ wantGrid?: boolean;
28+ sources: MatlabSources;
29+}
30+
31+export interface RunPoint {
32+ n: number;
33+ solveSeconds: number;
34+ solveSecondsAll: number[];
35+ coldSeconds: number;
36+ relMax: number;
37+ relL2: number;
38+ uEval: Float64Array;
39+ uGrid: Float64Array | null;
40+}
41+
42+function numLiteral(x: number): string {
43+ if (!Number.isFinite(x)) throw new Error(`bad numeric literal: ${x}`);
44+ return String(x);
45+}
46+
47+export function runPoint(req: RunPointRequest): RunPoint {
48+ const { instance, n } = req;
49+ const repeats = req.repeats ?? 3;
50+ const wantGrid = req.wantGrid ?? false;
51+ const main = [
52+ "% generated by the fastandaccurate harness",
53+ `prob = build_problem(${numLiteral(instance.a)}, ${numLiteral(instance.k)}, ${numLiteral(instance.d)}, ${wantGrid ? 1 : 0});`,
54+ `tic; out = solver(prob, ${n}); res_cold = toc;`,
55+ `res_times = zeros(${repeats}, 1);`,
56+ `for irep = 1:${repeats}`,
57+ ` tic; out = solver(prob, ${n}); res_times(irep) = toc;`,
58+ "end",
59+ "res_ueval = out.uEval;",
60+ "res_ugrid = out.uGrid;",
61+ "",
62+ ].join("\n");
63+ const { vars } = runNumblScript(
64+ main,
65+ {
66+ "build_problem.m": req.sources.buildProblem,
67+ "laplace2d_bdata.m": req.sources.bdata,
68+ "solver.m": req.sources.solver,
69+ },
70+ ["res_cold", "res_times", "res_ueval", "res_ugrid"]
71+ );
72+ const times = Array.from(vars.res_times);
73+ const sorted = [...times].sort((x, y) => x - y);
74+ const median =
75+ sorted.length % 2 === 1
76+ ? sorted[(sorted.length - 1) / 2]
77+ : (sorted[sorted.length / 2 - 1] + sorted[sorted.length / 2]) / 2;
78+ const { relMax, relL2 } = evalErrors(instance, vars.res_ueval);
79+ return {
80+ n,
81+ solveSeconds: median,
82+ solveSecondsAll: times,
83+ coldSeconds: vars.res_cold[0],
84+ relMax,
85+ relL2,
86+ uEval: vars.res_ueval,
87+ uGrid: wantGrid ? vars.res_ugrid : null,
88+ };
89+}
src/harness/sweep.tsadded+34−0View file
@@ -0,0 +1,34 @@
1+// A work-precision sweep: run one solver across its resolution list on
2+// one instance, reporting each point as it lands.
3+
4+import { runPoint, type MatlabSources, type RunPoint } from "./runner";
5+import type { Laplace2dInstance } from "../problems/laplace2d/spec";
6+import type { SolverManifest } from "../solvers";
7+
8+export interface SweepOptions {
9+ instance: Laplace2dInstance;
10+ solver: SolverManifest;
11+ sources: MatlabSources;
12+ repeats?: number;
13+ /** Restrict the sweep to n values <= this (for quick runs). */
14+ maxN?: number;
15+ onPoint?: (point: RunPoint, index: number, total: number) => void;
16+}
17+
18+export function runSweep(opts: SweepOptions): RunPoint[] {
19+ const ns = opts.solver.sweepN.filter(
20+ (n) => opts.maxN === undefined || n <= opts.maxN
21+ );
22+ const points: RunPoint[] = [];
23+ ns.forEach((n, i) => {
24+ const p = runPoint({
25+ instance: opts.instance,
26+ n,
27+ repeats: opts.repeats ?? 3,
28+ sources: opts.sources,
29+ });
30+ points.push(p);
31+ opts.onPoint?.(p, i, ns.length);
32+ });
33+ return points;
34+}
src/problems/laplace2d/exact.tsadded+125−0View file
@@ -0,0 +1,125 @@
1+// Exact solution, geometry, and sampling rules for laplace-dirichlet-2d,
2+// implemented independently of the MATLAB side (build_problem.m). The two
3+// implementations of these simple formulas check each other: a solver can
4+// only reach high accuracy if both agree.
5+
6+import type { Laplace2dInstance } from "./spec";
7+
8+export interface Source {
9+ x: number;
10+ y: number;
11+ c: number;
12+}
13+
14+/** Boundary radius r(t) = 1 + a cos(k t). */
15+export function boundaryR(inst: Laplace2dInstance, t: number): number {
16+ return 1 + inst.a * Math.cos(inst.k * t);
17+}
18+
19+/** Boundary point at parameter t. */
20+export function boundaryPoint(inst: Laplace2dInstance, t: number) {
21+ const r = boundaryR(inst, t);
22+ return { x: r * Math.cos(t), y: r * Math.sin(t) };
23+}
24+
25+/** The three exact-solution sources: boundary points at phi_j pushed a
26+ * distance d along the outward normal, strengths [1.0, -0.6, 0.8]. */
27+export function sources(inst: Laplace2dInstance): Source[] {
28+ const { a, k, d } = inst;
29+ const strengths = [1.0, -0.6, 0.8];
30+ return strengths.map((c, j) => {
31+ const phi = (2 * Math.PI * j) / 3 + 0.4;
32+ const r = 1 + a * Math.cos(k * phi);
33+ const dr = -a * k * Math.sin(k * phi);
34+ const bx = r * Math.cos(phi);
35+ const by = r * Math.sin(phi);
36+ const dx = dr * Math.cos(phi) - r * Math.sin(phi);
37+ const dy = dr * Math.sin(phi) + r * Math.cos(phi);
38+ const sp = Math.hypot(dx, dy);
39+ return { x: bx + (d * dy) / sp, y: by - (d * dx) / sp, c };
40+ });
41+}
42+
43+/** Exact solution u(x, y) = sum_j c_j log|x - s_j|. */
44+export function exactU(inst: Laplace2dInstance, x: number, y: number): number {
45+ let u = 0;
46+ for (const s of sources(inst)) {
47+ u += s.c * 0.5 * Math.log((x - s.x) ** 2 + (y - s.y) ** 2);
48+ }
49+ return u;
50+}
51+
52+/** The 65 evaluation points: 16 rays, radial fractions
53+ * [0.25, 0.5, 0.75, 0.9], plus the origin. Order matches
54+ * build_problem.m: radius outer, angle inner, origin last. */
55+export function evalPoints(inst: Laplace2dInstance): { x: number; y: number }[] {
56+ const rho = [0.25, 0.5, 0.75, 0.9];
57+ const pts: { x: number; y: number }[] = [];
58+ for (const r of rho) {
59+ for (let j = 0; j < 16; j++) {
60+ const th = (2 * Math.PI * j) / 16 + 0.13;
61+ const rr = r * boundaryR(inst, th);
62+ pts.push({ x: rr * Math.cos(th), y: rr * Math.sin(th) });
63+ }
64+ }
65+ pts.push({ x: 0, y: 0 });
66+ return pts;
67+}
68+
69+/** Exact solution at the evaluation points. */
70+export function exactAtEvalPoints(inst: Laplace2dInstance): Float64Array {
71+ const pts = evalPoints(inst);
72+ const u = new Float64Array(pts.length);
73+ pts.forEach((p, i) => (u[i] = exactU(inst, p.x, p.y)));
74+ return u;
75+}
76+
77+/** The visualization grid: ngrid x ngrid points over [-R, R]^2 with
78+ * R = 1.05 (1 + |a|). Flat index p = ix * ngrid + iy with x = xs[ix],
79+ * y = xs[iy] (y varies fastest), matching build_problem.m's meshgrid
80+ * column order. */
81+export const VIZ_NGRID = 200;
82+
83+export function vizGrid(inst: Laplace2dInstance) {
84+ const R = 1.05 * (1 + Math.abs(inst.a));
85+ const xs = new Float64Array(VIZ_NGRID);
86+ for (let i = 0; i < VIZ_NGRID; i++) {
87+ xs[i] = -R + (2 * R * i) / (VIZ_NGRID - 1);
88+ }
89+ return { R, ngrid: VIZ_NGRID, xs };
90+}
91+
92+/** Whether (x, y) is inside the domain (used only for display masking,
93+ * never for scoring, so float tie-breaks at the boundary are harmless). */
94+export function insideDomain(inst: Laplace2dInstance, x: number, y: number): boolean {
95+ const rr = Math.hypot(x, y);
96+ const th = Math.atan2(y, x);
97+ return rr < boundaryR(inst, th);
98+}
99+
100+/** Relative errors of numeric values against the exact solution at the
101+ * evaluation points: max and L2, both relative to the exact values. */
102+export function evalErrors(
103+ inst: Laplace2dInstance,
104+ uNum: ArrayLike<number>
105+): { relMax: number; relL2: number } {
106+ const uEx = exactAtEvalPoints(inst);
107+ if (uNum.length !== uEx.length) {
108+ throw new Error(`expected ${uEx.length} values, got ${uNum.length}`);
109+ }
110+ let maxDiff = 0;
111+ let maxEx = 0;
112+ let sumDiff2 = 0;
113+ let sumEx2 = 0;
114+ for (let i = 0; i < uEx.length; i++) {
115+ const diff = Math.abs(uNum[i] - uEx[i]);
116+ maxDiff = Math.max(maxDiff, diff);
117+ maxEx = Math.max(maxEx, Math.abs(uEx[i]));
118+ sumDiff2 += diff * diff;
119+ sumEx2 += uEx[i] * uEx[i];
120+ }
121+ return {
122+ relMax: maxDiff / maxEx,
123+ relL2: Math.sqrt(sumDiff2 / sumEx2),
124+ };
125+}
src/problems/laplace2d/matlab/build_problem.madded+73−0View file
@@ -0,0 +1,73 @@
1+function prob = build_problem(a, k, d, wantGrid)
2+% BUILD_PROBLEM Assemble the problem struct for laplace-dirichlet-2d.
3+%
4+% The domain is the star-shaped region bounded by
5+% x(t) = r(t) [cos t; sin t], r(t) = 1 + a cos(k t), t in [0, 2 pi).
6+% The exact solution is u(x) = sum_j c_j log|x - s_j|, with three point
7+% sources s_j outside the domain: s_j is the boundary point at parameter
8+% phi_j = 2 pi (j-1)/3 + 0.4 pushed a distance d along the outward unit
9+% normal, with strengths c = [1.0; -0.6; 0.8].
10+%
11+% The solver receives only the curve (with derivatives), the Dirichlet
12+% data g as a function of the boundary parameter, and the points where
13+% the solution is requested. The sources exist here only to manufacture
14+% the data; a submitted solver must not use knowledge of them.
15+%
16+% Fields of prob:
17+% curve @(t) -> [x y] boundary point, t column vector
18+% curveD @(t) -> [x' y'] first derivative
19+% curveDD @(t) -> [x'' y''] second derivative
20+% g @(t) -> g Dirichlet data at boundary parameter t
21+% evalXY 65 x 2 points where uEval is required
22+% vizXY m x 2 grid points where uGrid is requested
23+% (m = 0 when no visualization is wanted)
24+
25+phi = 2*pi*[0; 1; 2]/3 + 0.4;
26+c = [1.0; -0.6; 0.8];
27+rphi = 1 + a*cos(k*phi);
28+bx = rphi.*cos(phi);
29+by = rphi.*sin(phi);
30+dxb = -a*k*sin(k*phi).*cos(phi) - rphi.*sin(phi);
31+dyb = -a*k*sin(k*phi).*sin(phi) + rphi.*cos(phi);
32+sp = sqrt(dxb.^2 + dyb.^2);
33+sx = bx + d*(dyb./sp);
34+sy = by - d*(dxb./sp);
35+
36+prob = struct();
37+prob.curve = @(t) [(1 + a*cos(k*t)).*cos(t), (1 + a*cos(k*t)).*sin(t)];
38+prob.curveD = @(t) [-a*k*sin(k*t).*cos(t) - (1 + a*cos(k*t)).*sin(t), ...
39+ -a*k*sin(k*t).*sin(t) + (1 + a*cos(k*t)).*cos(t)];
40+prob.curveDD = @(t) [(-a*k*k*cos(k*t) - 1 - a*cos(k*t)).*cos(t) + 2*a*k*sin(k*t).*sin(t), ...
41+ (-a*k*k*cos(k*t) - 1 - a*cos(k*t)).*sin(t) - 2*a*k*sin(k*t).*cos(t)];
42+prob.g = @(t) laplace2d_bdata(t, a, k, sx, sy, c);
43+
44+% Evaluation points: 16 rays, 4 radial fractions, plus the origin.
45+% The rule must match evalPoints() in src/problems/laplace2d/exact.ts.
46+rho = [0.25; 0.5; 0.75; 0.9];
47+th = 2*pi*(0:15)'/16 + 0.13;
48+pts = zeros(numel(rho)*numel(th) + 1, 2);
49+idx = 1;
50+for i = 1:numel(rho)
51+ for j = 1:numel(th)
52+ rr = rho(i)*(1 + a*cos(k*th(j)));
53+ pts(idx, 1) = rr*cos(th(j));
54+ pts(idx, 2) = rr*sin(th(j));
55+ idx = idx + 1;
56+ end
57+end
58+prob.evalXY = pts;
59+
60+% Visualization grid: ngrid x ngrid points over the bounding square,
61+% listed with y varying fastest (MATLAB column order). Points outside
62+% the domain are included; the viewer masks them.
63+if wantGrid
64+ ngrid = 200;
65+ R = 1.05*(1 + abs(a));
66+ xs = linspace(-R, R, ngrid);
67+ [X, Y] = meshgrid(xs, xs);
68+ prob.vizXY = [X(:), Y(:)];
69+else
70+ prob.vizXY = zeros(0, 2);
71+end
72+
73+end
src/problems/laplace2d/matlab/laplace2d_bdata.madded+13−0View file
@@ -0,0 +1,13 @@
1+function g = laplace2d_bdata(t, a, k, sx, sy, c)
2+% LAPLACE2D_BDATA Dirichlet data for laplace-dirichlet-2d.
3+% Evaluates u(x(t)) = sum_j c_j log|x(t) - s_j| at boundary parameters t.
4+% Called through the prob.g handle built in build_problem; solvers see
5+% only that handle.
6+r = 1 + a*cos(k*t);
7+x = r.*cos(t);
8+y = r.*sin(t);
9+g = zeros(size(t));
10+for j = 1:numel(c)
11+ g = g + c(j)*0.5*log((x - sx(j)).^2 + (y - sy(j)).^2);
12+end
13+end
src/problems/laplace2d/spec.tsadded+83−0View file
@@ -0,0 +1,83 @@
1+// Problem: laplace-dirichlet-2d.
2+// The canonical statement lives in docs/problems/laplace-dirichlet-2d.md.
3+// This module defines the official instances and the canonical spec object
4+// that identifies a (problem, instance) pair in results and, later, in
5+// cache keys.
6+
7+export const PROBLEM_ID = "laplace-dirichlet-2d";
8+export const PROBLEM_VERSION = 1;
9+
10+export interface Laplace2dInstance {
11+ /** Short stable identifier used in results and URLs. */
12+ id: string;
13+ label: string;
14+ /** Boundary r(t) = 1 + a cos(k t). */
15+ a: number;
16+ k: number;
17+ /** Distance of the exact solution's sources beyond the boundary. */
18+ d: number;
19+ description: string;
20+}
21+
22+export const INSTANCES: Laplace2dInstance[] = [
23+ {
24+ id: "disk-easy",
25+ label: "Disk, distant sources",
26+ a: 0,
27+ k: 0,
28+ d: 0.5,
29+ description:
30+ "The unit disk with sources half a radius beyond the boundary. " +
31+ "Every reasonable method should reach high accuracy quickly.",
32+ },
33+ {
34+ id: "star-medium",
35+ label: "3-lobe star, moderate sources",
36+ a: 0.2,
37+ k: 3,
38+ d: 0.4,
39+ description:
40+ "A gently star-shaped domain; the data continues comfortably past " +
41+ "the boundary, so geometric convergence is attainable but the " +
42+ "geometry is no longer trivial.",
43+ },
44+ {
45+ id: "star-hard",
46+ label: "5-lobe star, close sources",
47+ a: 0.3,
48+ k: 5,
49+ d: 0.08,
50+ description:
51+ "A wavier domain with sources only 0.08 beyond the boundary. The " +
52+ "data barely continues past the boundary, which defeats methods " +
53+ "whose representation assumes it does.",
54+ },
55+];
56+
57+export function getInstance(id: string): Laplace2dInstance {
58+ const inst = INSTANCES.find((i) => i.id === id);
59+ if (!inst) throw new Error(`Unknown instance: ${id}`);
60+ return inst;
61+}
62+
63+/**
64+ * The canonical spec object for an instance. Serialized with sorted keys,
65+ * this string identifies the instance exactly (results carry it, and a
66+ * future artifact cache hashes it).
67+ */
68+export function canonicalSpec(inst: Laplace2dInstance) {
69+ return {
70+ a: inst.a,
71+ d: inst.d,
72+ instance: inst.id,
73+ k: inst.k,
74+ problem: PROBLEM_ID,
75+ problemVersion: PROBLEM_VERSION,
76+ };
77+}
78+
79+export function canonicalSpecJson(inst: Laplace2dInstance): string {
80+ const spec = canonicalSpec(inst);
81+ const keys = Object.keys(spec).sort();
82+ return JSON.stringify(spec, keys);
83+}
src/solvers/index.tsadded+54−0View file
@@ -0,0 +1,54 @@
1+// The in-browser solvers shipped with the site. Each solver is a MATLAB
2+// function file (solver.m) implementing the interface documented in
3+// docs/problems/laplace-dirichlet-2d.md. A submitted solver adds a
4+// directory here plus a manifest entry.
5+
6+export interface SolverManifest {
7+ /** Short stable identifier used in results and URLs. */
8+ id: string;
9+ name: string;
10+ /** One-paragraph description of the method for the problem page. */
11+ description: string;
12+ /** Semantic version of the solver code; bump on any change that could
13+ * alter results. */
14+ version: string;
15+ backend: "cpu" | "gpu";
16+ /** The resolution values a standard work-precision sweep runs. */
17+ sweepN: number[];
18+}
19+
20+export const SOLVERS: SolverManifest[] = [
21+ {
22+ id: "mfs",
23+ name: "Method of fundamental solutions",
24+ description:
25+ "Represents the solution as n logarithmic point charges on a curve " +
26+ "a fixed distance 0.3 outside the boundary, with strengths found by " +
27+ "collocation at n boundary points. Converges geometrically when the " +
28+ "data continues harmonically past the charge curve; stagnates when " +
29+ "it does not (the star-hard instance). Ill-conditioning caps the " +
30+ "attainable accuracy near 1e-10 in exchange for very small n.",
31+ version: "1.0.0",
32+ backend: "cpu",
33+ sweepN: [8, 12, 16, 24, 32, 48, 64, 96, 128, 192, 256],
34+ },
35+ {
36+ id: "nystrom-dlp",
37+ name: "Nystrom double-layer BIE",
38+ description:
39+ "Second-kind boundary integral equation for the double-layer " +
40+ "density, discretized with the periodic trapezoid rule at n " +
41+ "boundary nodes. Converges geometrically for any smooth data, at a " +
42+ "rate set by how far the data continues analytically, so the hard " +
43+ "instance costs more nodes rather than a lost method assumption.",
44+ version: "1.0.0",
45+ backend: "cpu",
46+ sweepN: [16, 24, 32, 48, 64, 96, 128, 192, 256, 384, 512, 768],
47+ },
48+];
49+
50+export function getSolver(id: string): SolverManifest {
51+ const s = SOLVERS.find((x) => x.id === id);
52+ if (!s) throw new Error(`Unknown solver: ${id}`);
53+ return s;
54+}
src/solvers/mfs/solver.madded+60−0View file
@@ -0,0 +1,60 @@
1+function out = solver(prob, n)
2+% Method of fundamental solutions (MFS) for the interior Dirichlet
3+% Laplace problem.
4+%
5+% The solution is represented as a sum of n logarithmic point charges
6+% placed on a fictitious curve outside the domain: each charge sits at
7+% the boundary point of parameter t_j pushed a fixed distance delta
8+% along the outward normal. Charge strengths are found by collocation
9+% of the Dirichlet data at the n boundary points t_j (a square dense
10+% system). Convergence is geometric when the solution continues
11+% harmonically past the charge curve; when the data has singularities
12+% closer to the boundary than delta, the method stagnates. The system
13+% grows exponentially ill-conditioned with n, which caps the attainable
14+% accuracy near 1e-10.
15+%
16+% n : number of charges = number of collocation points.
17+
18+delta = 0.3;
19+
20+t = 2*pi*(0:n-1)'/n;
21+xb = prob.curve(t);
22+dxb = prob.curveD(t);
23+sp = sqrt(dxb(:,1).^2 + dxb(:,2).^2);
24+nx = dxb(:,2)./sp;
25+ny = -dxb(:,1)./sp;
26+qx = xb(:,1) + delta*nx;
27+qy = xb(:,2) + delta*ny;
28+
29+A = kernel_matrix(xb(:,1), xb(:,2), qx, qy);
30+coef = A \ prob.g(t);
31+
32+out = struct();
33+out.uEval = apply_potential(prob.evalXY, qx, qy, coef);
34+if size(prob.vizXY, 1) > 0
35+ out.uGrid = apply_potential(prob.vizXY, qx, qy, coef);
36+else
37+ out.uGrid = zeros(0, 1);
38+end
39+
40+end
41+
42+function A = kernel_matrix(px, py, qx, qy)
43+% A(i,j) = log|p_i - q_j|
44+m = numel(px);
45+nq = numel(qx);
46+dx = repmat(px, 1, nq) - repmat(qx', m, 1);
47+dy = repmat(py, 1, nq) - repmat(qy', m, 1);
48+A = 0.5*log(dx.^2 + dy.^2);
49+end
50+
51+function u = apply_potential(XY, qx, qy, coef)
52+% Evaluate the charge sum at the rows of XY, in blocks to bound memory.
53+m = size(XY, 1);
54+u = zeros(m, 1);
55+B = 4000;
56+for i0 = 1:B:m
57+ i1 = min(i0 + B - 1, m);
58+ u(i0:i1) = kernel_matrix(XY(i0:i1, 1), XY(i0:i1, 2), qx, qy)*coef;
59+end
60+end
src/solvers/nystrom-dlp/solver.madded+70−0View file
@@ -0,0 +1,70 @@
1+function out = solver(prob, n)
2+% Nystrom discretization of the double-layer boundary integral equation
3+% for the interior Dirichlet Laplace problem.
4+%
5+% The solution is represented as a double-layer potential
6+% u(x) = (1/2pi) int_Gamma sigma(y) (x - y).n(y) / |x - y|^2 ds(y)
7+% whose interior boundary limit gives the second-kind equation
8+% (W - I/2) sigma = g.
9+% W is discretized with the periodic trapezoid rule at n equispaced
10+% parameter nodes; the kernel is smooth on a smooth curve, with the
11+% diagonal limit -kappa(t) |x'(t)| / (4 pi). Convergence is geometric,
12+% at a rate set by how far the data g continues analytically. Accuracy
13+% of the evaluated potential degrades for targets very close to the
14+% boundary (the close-evaluation problem); the evaluation points of
15+% this problem stay a modest distance inside.
16+%
17+% n : number of boundary quadrature nodes.
18+
19+h = 2*pi/n;
20+t = h*(0:n-1)';
21+xb = prob.curve(t);
22+dxb = prob.curveD(t);
23+ddxb = prob.curveDD(t);
24+sp = sqrt(dxb(:,1).^2 + dxb(:,2).^2);
25+nqx = dxb(:,2)./sp;
26+nqy = -dxb(:,1)./sp;
27+
28+% M(i,j) = (1/2pi) (x_i - x_j).n(x_j) / |x_i - x_j|^2 * |x'(t_j)|
29+dx = repmat(xb(:,1), 1, n) - repmat(xb(:,1)', n, 1);
30+dy = repmat(xb(:,2), 1, n) - repmat(xb(:,2)', n, 1);
31+r2 = dx.^2 + dy.^2;
32+num = dx.*repmat(nqx', n, 1) + dy.*repmat(nqy', n, 1);
33+M = (num./r2).*repmat(sp', n, 1)/(2*pi);
34+
35+% Diagonal limit: -kappa/2 * |x'| / (2 pi), kappa the signed curvature.
36+kap = (dxb(:,1).*ddxb(:,2) - dxb(:,2).*ddxb(:,1))./sp.^3;
37+md = -(kap/2).*sp/(2*pi);
38+for i = 1:n
39+ M(i, i) = md(i);
40+end
41+
42+sigma = (h*M - 0.5*eye(n)) \ prob.g(t);
43+
44+out = struct();
45+out.uEval = dlp_eval(prob.evalXY, xb, nqx, nqy, h*sp.*sigma);
46+if size(prob.vizXY, 1) > 0
47+ out.uGrid = dlp_eval(prob.vizXY, xb, nqx, nqy, h*sp.*sigma);
48+else
49+ out.uGrid = zeros(0, 1);
50+end
51+
52+end
53+
54+function u = dlp_eval(XY, xb, nqx, nqy, w)
55+% Evaluate the double-layer potential with combined weights w at the
56+% rows of XY, in blocks to bound memory.
57+m = size(XY, 1);
58+nb = size(xb, 1);
59+u = zeros(m, 1);
60+B = 4000;
61+for i0 = 1:B:m
62+ i1 = min(i0 + B - 1, m);
63+ mm = i1 - i0 + 1;
64+ dx = repmat(XY(i0:i1, 1), 1, nb) - repmat(xb(:,1)', mm, 1);
65+ dy = repmat(XY(i0:i1, 2), 1, nb) - repmat(xb(:,2)', mm, 1);
66+ r2 = dx.^2 + dy.^2;
67+ K = (dx.*repmat(nqx', mm, 1) + dy.*repmat(nqy', mm, 1))./r2/(2*pi);
68+ u(i0:i1) = K*w;
69+end
70+end
test/solver-test.tsadded+95−0View file
@@ -0,0 +1,95 @@
1+// Convergence test: run both solvers on the official instances through
2+// numbl in node and check that errors behave as the theory says they
3+// should. Run with: npx tsx test/solver-test.ts
4+
5+import { readFileSync } from "fs";
6+import { fileURLToPath } from "url";
7+import { dirname, join } from "path";
8+import { INSTANCES, getInstance } from "../src/problems/laplace2d/spec";
9+import { SOLVERS } from "../src/solvers";
10+import { runPoint, type MatlabSources } from "../src/harness/runner";
11+import { runSweep } from "../src/harness/sweep";
12+
13+const root = join(dirname(fileURLToPath(import.meta.url)), "..");
14+const read = (p: string) => readFileSync(join(root, p), "utf-8");
15+
16+const base = {
17+ buildProblem: read("src/problems/laplace2d/matlab/build_problem.m"),
18+ bdata: read("src/problems/laplace2d/matlab/laplace2d_bdata.m"),
19+};
20+const solverSources: Record<string, MatlabSources> = {
21+ mfs: { ...base, solver: read("src/solvers/mfs/solver.m") },
22+ "nystrom-dlp": { ...base, solver: read("src/solvers/nystrom-dlp/solver.m") },
23+};
24+
25+// Best relMax each solver must reach over its full sweep. On star-hard,
26+// MFS is additionally required NOT to do well: its charge curve lies
27+// beyond the data's singularities there, and if it suddenly reached high
28+// accuracy the instance would no longer be testing what the spec says.
29+const mustReach: Record<string, Record<string, number>> = {
30+ mfs: { "disk-easy": 1e-12, "star-medium": 1e-12, "star-hard": 1e-2 },
31+ "nystrom-dlp": { "disk-easy": 1e-10, "star-medium": 1e-10, "star-hard": 1e-8 },
32+};
33+const mustNotReach: Record<string, Record<string, number>> = {
34+ mfs: { "star-hard": 1e-8 },
35+};
36+
37+let failures = 0;
38+
39+for (const inst of INSTANCES) {
40+ for (const solver of SOLVERS) {
41+ console.log(`\n== ${inst.id} / ${solver.id}`);
42+ console.log(" n relMax relL2 solve(s) cold(s)");
43+ let best = Infinity;
44+ runSweep({
45+ instance: inst,
46+ solver,
47+ sources: solverSources[solver.id],
48+ repeats: 1,
49+ onPoint: (p) => {
50+ best = Math.min(best, p.relMax);
51+ console.log(
52+ ` ${String(p.n).padStart(4)} ${p.relMax.toExponential(3)} ` +
53+ `${p.relL2.toExponential(3)} ${p.solveSeconds.toFixed(4)} ${p.coldSeconds.toFixed(4)}`
54+ );
55+ },
56+ });
57+ const reach = mustReach[solver.id][inst.id];
58+ const notReach = mustNotReach[solver.id]?.[inst.id];
59+ if (best > reach) {
60+ console.log(` FAIL: best relMax ${best.toExponential(2)} > ${reach}`);
61+ failures++;
62+ } else if (notReach !== undefined && best < notReach) {
63+ console.log(
64+ ` FAIL: best relMax ${best.toExponential(2)} < ${notReach} ` +
65+ "(instance no longer defeats this method)"
66+ );
67+ failures++;
68+ } else {
69+ console.log(` ok (best relMax ${best.toExponential(2)})`);
70+ }
71+ }
72+}
73+
74+// The grid path: one run with wantGrid on the medium instance.
75+{
76+ const p = runPoint({
77+ instance: getInstance("star-medium"),
78+ n: 64,
79+ repeats: 1,
80+ wantGrid: true,
81+ sources: solverSources["nystrom-dlp"],
82+ });
83+ if (!p.uGrid || p.uGrid.length !== 200 * 200) {
84+ console.log(`\nFAIL: grid has ${p.uGrid?.length ?? 0} values, expected 40000`);
85+ failures++;
86+ } else {
87+ console.log(`\ngrid ok (${p.uGrid.length} values)`);
88+ }
89+}
90+
91+if (failures > 0) {
92+ console.error(`\n${failures} failure(s)`);
93+ process.exit(1);
94+}
95+console.log("\nall checks passed");
tsconfig.jsonadded+17−0View file
@@ -0,0 +1,17 @@
1+{
2+ "compilerOptions": {
3+ "target": "ES2022",
4+ "lib": ["ES2022", "DOM", "DOM.Iterable", "WebWorker"],
5+ "module": "ESNext",
6+ "moduleResolution": "bundler",
7+ "jsx": "react-jsx",
8+ "strict": true,
9+ "noEmit": true,
10+ "skipLibCheck": true,
11+ "isolatedModules": true,
12+ "noUnusedLocals": true,
13+ "noUnusedParameters": true,
14+ "types": ["vite/client", "node"]
15+ },
16+ "include": ["src", "test", "scripts", "vite.config.ts"]
17+}
vite.cli.config.tsadded+32−0View file
@@ -0,0 +1,32 @@
1+/**
2+ * The command line's bundle: the same harness source as the page, built
3+ * for node. numbl is pure TypeScript/JavaScript, so everything bundles
4+ * and the published package has no dependencies at all. The MATLAB
5+ * sources are not bundled; pack-cli.mjs copies them into the tarball and
6+ * the CLI reads them from disk beside itself, so the installed package
7+ * shows the exact solver code it runs.
8+ */
9+import { defineConfig, mergeConfig } from "vite";
10+import base from "./vite.config";
11+
12+export default mergeConfig(
13+ base,
14+ defineConfig({
15+ ssr: {
16+ // Bundle everything (numbl included); the published tarball has no
17+ // dependencies.
18+ noExternal: true,
19+ },
20+ build: {
21+ ssr: "src/cli/main.ts",
22+ outDir: "dist-cli",
23+ target: "node20",
24+ emptyOutDir: true,
25+ minify: false,
26+ rollupOptions: {
27+ input: "src/cli/main.ts",
28+ output: { entryFileNames: "main.js" },
29+ },
30+ },
31+ })
32+);
vite.config.tsadded+26−0View file
@@ -0,0 +1,26 @@
1+import { defineConfig } from "vite";
2+import react from "@vitejs/plugin-react";
3+import { readFileSync } from "fs";
4+import { fileURLToPath } from "url";
5+import { dirname, join } from "path";
6+
7+const root = dirname(fileURLToPath(import.meta.url));
8+const numblVersion = (
9+ JSON.parse(
10+ readFileSync(join(root, "node_modules", "numbl", "package.json"), "utf-8")
11+ ) as { version: string }
12+).version;
13+
14+export default defineConfig({
15+ base: "./",
16+ plugins: [react()],
17+ define: {
18+ __NUMBL_VERSION__: JSON.stringify(numblVersion),
19+ __BUILD_ID__: JSON.stringify(
20+ process.env.GITHUB_SHA?.slice(0, 7) ?? "dev"
21+ ),
22+ },
23+ worker: {
24+ format: "es",
25+ },
26+});