Add chunkie-dlp solver: mip package support in the harness, curl-backed file I/O for node
17 changed files+247−212
docs/problems/laplace-dirichlet-2d.mdmodified+16−50View file
@@ -21,52 +21,28 @@ set of evaluation points.
2121
2222 ## The exact solution
2323
24-The data g is manufactured from an exact harmonic function u*, and
25-g = u* restricted to ∂Ω. Since u* is harmonic in Ω, it is the unique
26-solution, and errors are measured against it directly rather than
27-against a reference computation. Two data families are used.
28-
29-**log-sources.** A sum of three logarithmic point sources placed outside
30-the domain:
24+The data g is manufactured from an exact harmonic function, a sum of
25+three logarithmic point sources placed outside the domain:
3126
3227 u*(x) = Σ_j c_j log |x − s_j|, c = (1.0, −0.6, 0.8).
3328
3429 Source s_j is the boundary point at parameter φ_j = 2π(j−1)/3 + 0.4
35-pushed a distance d along the outward unit normal.
36-
37-**branch-point.** The real part of a complex square root:
38-
39- u*(z) = Re √(w), w = −(z − z₀) e^{−iθ₀},
40-
41-where z₀ is the boundary point at parameter 0.4 pushed a distance d
42-along the outward unit normal and θ₀ is its polar angle. The principal
43-branch cut w ≤ 0 maps to the radial ray from z₀ away from the origin,
44-which cannot meet a domain star-shaped about the origin, so u* is
45-harmonic on the closed domain; its nearest singularity is the branch
46-point z₀ itself. In real arithmetic, Re √w = √((|w| + Re w)/2).
47-
48-In both families the distance d controls difficulty: u* continues
49-harmonically only up to its singularities, so the smaller d, the shorter
50-the distance the data continues past the boundary, and methods whose
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
5136 representations assume a generous continuation lose it. A solver must
52-not use knowledge of the sources or branch point; they exist only to
53-manufacture g.
37+not use knowledge of the sources; they exist only to manufacture g.
5438
5539 ## Official instances
5640
57-| id | a | k | d | family | character |
58-|---|---|---|---|---|---|
59-| disk-easy | 0 | 0 | 0.5 | log-sources | the unit disk, distant sources |
60-| star-medium | 0.2 | 3 | 0.4 | log-sources | mild geometry, comfortable continuation |
61-| star-hard | 0.3 | 5 | 0.08 | log-sources | wavy boundary, data barely continues |
62-| star-branch | 0.2 | 3 | 0.4 | branch-point | star-medium's geometry with the other data family |
63-
64-star-branch repeats star-medium's geometry and singularity distance with
65-the branch-point family. The pair separates the data family from its
66-singularity structure: a method whose accuracy depended on a coincidence
67-between its representation and the data-generating family would score
68-differently on the two, while a method governed by the continuation
69-distance scores alike.
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 |
7046
7147 Committed results exist only at these instances, so every solver is
7248 compared on identical inputs. The parameters may be varied freely in the
@@ -140,15 +116,5 @@ masks them; grid values are never scored.
140116 The exact solution is smooth and free of boundary singularities, so this
141117 problem does not test corner handling, nonsmooth data, or interior
142118 sources (a nonzero right-hand side would exclude plain boundary-integral
143-methods; that belongs to a different problem). In particular, every
144-instance's solution is real-analytic in a neighborhood of the closed
145-domain, which favors methods that exploit analyticity (fundamental
146-solutions, boundary integral equations, spectral collocation) over
147-low-order discretizations; this is a property of manufactured smooth
148-benchmarks in general, not of the particular families chosen. The
149-star-branch instance exists because most instances manufacture the data
150-from exterior log charges, the same function class the MFS
151-representation draws on; comparing it with star-medium checks that
152-scores follow the continuation distance rather than the family. The
153-domain family is star-shaped by construction, which some methods can
154-exploit.
119+methods; that belongs to a different problem). The domain family is
120+star-shaped by construction, which some methods can exploit.
package-lock.jsonmodified+1−0View file
@@ -8,6 +8,7 @@
88 "name": "fastandaccurate",
99 "version": "0.1.0",
1010 "dependencies": {
11+ "fflate": "^0.8.3",
1112 "numbl": "^0.4.18",
1213 "react": "^19.2.7",
1314 "react-dom": "^19.2.7"
package.jsonmodified+1−0View file
@@ -12,6 +12,7 @@
1212 "check-app": "node scripts/check-app.mjs"
1313 },
1414 "dependencies": {
15+ "fflate": "^0.8.3",
1516 "numbl": "^0.4.18",
1617 "react": "^19.2.7",
1718 "react-dom": "^19.2.7"
src/app/components/DomainView.tsxmodified+5−26View file
@@ -5,9 +5,8 @@ import { useEffect, useRef } from "react";
55 import type { Laplace2dInstance } from "../../problems/laplace2d/spec";
66 import {
77 boundaryPoint,
8- branchPoint,
98 evalPoints,
10- singularities,
9+ sources,
1110 } from "../../problems/laplace2d/exact";
1211
1312 function token(name: string): string {
@@ -60,10 +59,10 @@ export function DomainView({ inst }: { inst: Laplace2dInstance }) {
6059 ctx.fill();
6160 }
6261
63- // the exact solution's singularities
62+ // sources
6463 ctx.strokeStyle = token("--series-2");
6564 ctx.lineWidth = 2;
66- for (const src of singularities(inst)) {
65+ for (const src of sources(inst)) {
6766 const cx = X(src.x);
6867 const cy = Y(src.y);
6968 ctx.beginPath();
@@ -73,22 +72,6 @@ export function DomainView({ inst }: { inst: Laplace2dInstance }) {
7372 ctx.lineTo(cx + 5, cy - 5);
7473 ctx.stroke();
7574 }
76-
77- // branch family: dashed segment along the cut (radially outward)
78- if (inst.family === "branch-point") {
79- const b = branchPoint(inst);
80- ctx.strokeStyle = token("--text-3");
81- ctx.lineWidth = 1.5;
82- ctx.setLineDash([4, 3]);
83- ctx.beginPath();
84- ctx.moveTo(X(b.x), Y(b.y));
85- ctx.lineTo(
86- X(b.x + 0.3 * Math.cos(b.theta0)),
87- Y(b.y + 0.3 * Math.sin(b.theta0))
88- );
89- ctx.stroke();
90- ctx.setLineDash([]);
91- }
9275 };
9376 draw();
9477 const mq = window.matchMedia("(prefers-color-scheme: dark)");
@@ -101,12 +84,8 @@ export function DomainView({ inst }: { inst: Laplace2dInstance }) {
10184 <canvas ref={canvasRef} />
10285 <figcaption className="field-caption" style={{ maxWidth: 360 }}>
10386 The domain, the evaluation points where solutions are scored
104- (dots), and the exact solution's{" "}
105- {inst.family === "branch-point"
106- ? "branch point a distance " + inst.d + " outside the boundary " +
107- "(cross), with its cut running radially outward (dashed)"
108- : "sources a distance " + inst.d + " outside the boundary (crosses)"}
109- .
87+ (dots), and the exact solution's sources a distance {inst.d} outside
88+ the boundary (crosses).
11089 </figcaption>
11190 </figure>
11291 );
src/app/matlabSources.tsmodified+0−1View file
@@ -17,7 +17,6 @@ export function matlabBase() {
1717 return {
1818 buildProblem: get("../problems/laplace2d/matlab/build_problem.m"),
1919 bdata: get("../problems/laplace2d/matlab/laplace2d_bdata.m"),
20- bdataBranch: get("../problems/laplace2d/matlab/laplace2d_bdata_branch.m"),
2120 };
2221 }
2322
src/app/pages/ProblemPage.tsxmodified+0−2View file
@@ -236,8 +236,6 @@ export function ProblemPage({ problemId }: { problemId: string }) {
236236 <td>{inst.k}</td>
237237 <th className="left">d</th>
238238 <td>{inst.d}</td>
239- <th className="left">data</th>
240- <td className="left">{inst.family}</td>
241239 </tr>
242240 </tbody>
243241 </table>
src/cli/main.tsmodified+4−1View file
@@ -15,6 +15,10 @@ import { dirname, join, resolve } from "path";
1515 import { fileURLToPath } from "url";
1616 import os from "os";
1717 import { INSTANCES, getInstance } from "../problems/laplace2d/spec";
18+import { setNumblFileIO } from "../harness/numblRun";
19+import { NodeFileIOAdapter } from "./nodeFileIO";
20+
21+setNumblFileIO((vfs) => new NodeFileIOAdapter(vfs));
1822 import { SOLVERS, getSolver, type SolverManifest } from "../solvers";
1923 import { runSweep } from "../harness/sweep";
2024 import type { MatlabSources } from "../harness/runner";
@@ -110,7 +114,6 @@ async function runCommand(flags: Record<string, string>) {
110114 const base = {
111115 buildProblem: readSrc("problems/laplace2d/matlab/build_problem.m"),
112116 bdata: readSrc("problems/laplace2d/matlab/laplace2d_bdata.m"),
113- bdataBranch: readSrc("problems/laplace2d/matlab/laplace2d_bdata_branch.m"),
114117 };
115118
116119 let solverList: { manifest: SolverManifest; sources: MatlabSources; source: string }[];
src/cli/nodeFileIO.tsadded+41−0View file
@@ -0,0 +1,41 @@
1+// The node substitute for numbl's synchronous-XHR file I/O: websave and
2+// webread shell out to curl (node has no synchronous fetch), with
3+// responses cached under ~/.cache/fastandaccurate keyed by URL, so
4+// repeated runs are offline. Used by the CLI and the node tests for
5+// solvers that mip-install packages.
6+
7+import { execFileSync } from "node:child_process";
8+import { createHash } from "node:crypto";
9+import { existsSync, mkdirSync, readFileSync } from "node:fs";
10+import { homedir } from "node:os";
11+import { join } from "node:path";
12+import { BrowserFileIOAdapter, type VirtualFileSystem } from "numbl";
13+
14+const cacheDir = join(homedir(), ".cache", "fastandaccurate");
15+
16+function curlCached(url: string): Buffer {
17+ mkdirSync(cacheDir, { recursive: true });
18+ const key = createHash("sha1").update(url).digest("hex").slice(0, 16);
19+ const cached = join(cacheDir, key);
20+ if (!existsSync(cached)) {
21+ execFileSync("curl", ["-sfL", "-o", cached, url], { stdio: "inherit" });
22+ }
23+ return readFileSync(cached);
24+}
25+
26+export class NodeFileIOAdapter extends BrowserFileIOAdapter {
27+ private nodeVfs: VirtualFileSystem;
28+ constructor(vfs: VirtualFileSystem) {
29+ super(vfs);
30+ this.nodeVfs = vfs;
31+ }
32+ override websave(url: string, filename: string): void {
33+ this.nodeVfs.writeFile(
34+ this.nodeVfs.normalizePath(filename),
35+ new Uint8Array(curlCached(url))
36+ );
37+ }
38+ override webread(url: string): string {
39+ return curlCached(url).toString("utf8");
40+ }
41+}
src/harness/numblRun.tsmodified+64−5View file
@@ -1,12 +1,68 @@
11 // 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.
2+// Works identically in a browser worker and in node. Plain solvers need
3+// no file system at all: the MATLAB side communicates through workspace
4+// variables, read back from result.variableValues.
5+//
6+// Solvers that begin with `mip load --install <pkg>` additionally need
7+// the mip package manager, which the wrapper bootstraps on first use: a
8+// persistent virtual file system holds /system, the mip core is fetched
9+// once and unzipped into it, and mip's own downloads go through the
10+// platform's websave. In a browser worker that is numbl's
11+// BrowserFileIOAdapter (synchronous XHR, GitHub release URLs routed
12+// through numbl's CORS proxy); the node CLI substitutes a curl-backed
13+// adapter via setNumblFileIO. Installed packages persist for the
14+// lifetime of the worker or process, so a sweep pays the download once.
515
6-import { executeCode } from "numbl";
16+import {
17+ executeCode,
18+ VirtualFileSystem,
19+ BrowserFileIOAdapter,
20+ BrowserSystemAdapter,
21+} from "numbl";
22+import { unzipSync } from "fflate";
723
824 const PROJ = "/fastandaccurate";
925
26+const MIP_MHL_URL =
27+ "https://github.com/mip-org/mip-core/releases/download/mip-numbl/mip-numbl-any.mhl";
28+const MIP_SYSTEM_PREFIX = "/system/mip/packages/gh/mip-org/core/mip/";
29+const MIP_SEARCH_PATH = MIP_SYSTEM_PREFIX + "mip";
30+
31+type FileIOFactory = (vfs: VirtualFileSystem) => BrowserFileIOAdapter;
32+
33+let makeFileIO: FileIOFactory = (vfs) => new BrowserFileIOAdapter(vfs);
34+
35+/** Substitute the platform's file I/O adapter (the node CLI installs a
36+ * curl-backed one; the browser default needs nothing). Must be called
37+ * before the first mip-using run. */
38+export function setNumblFileIO(factory: FileIOFactory) {
39+ makeFileIO = factory;
40+}
41+
42+let vfs: VirtualFileSystem | null = null;
43+let fileIO: BrowserFileIOAdapter | null = null;
44+let system: BrowserSystemAdapter | null = null;
45+let mipReady = false;
46+
47+function ensureMip() {
48+ if (mipReady && vfs && fileIO && system) return;
49+ vfs = new VirtualFileSystem();
50+ fileIO = makeFileIO(vfs);
51+ system = new BrowserSystemAdapter(vfs);
52+ const tmp = "/tmp/mip-core.mhl";
53+ fileIO.websave(MIP_MHL_URL, tmp);
54+ const entries = unzipSync(vfs.readFile(vfs.normalizePath(tmp)));
55+ for (const [name, content] of Object.entries(entries)) {
56+ if (name.endsWith("/")) continue;
57+ vfs.writeFile(MIP_SYSTEM_PREFIX + name, content);
58+ }
59+ mipReady = true;
60+}
61+
62+function usesMip(sources: string[]): boolean {
63+ return sources.some((s) => /^\s*mip\s+load\b/m.test(s));
64+}
65+
1066 export interface NumblRunResult {
1167 /** Console output of the run. */
1268 output: string;
@@ -30,6 +86,8 @@ export function runNumblScript(
3086 source,
3187 }));
3288 const outputs: string[] = [];
89+ const mip = usesMip([mainSource, ...Object.values(files)]);
90+ if (mip) ensureMip();
3391 const result = executeCode(
3492 mainSource,
3593 {
@@ -37,10 +95,11 @@ export function runNumblScript(
3795 displayResults: false,
3896 optimization: "1",
3997 implicitCwdPath: null,
98+ ...(mip && fileIO && system ? { fileIO, system } : {}),
4099 },
41100 workspaceFiles,
42101 `${PROJ}/main.m`,
43- [PROJ]
102+ mip ? [PROJ, MIP_SEARCH_PATH] : [PROJ]
44103 );
45104 const vars: Record<string, Float64Array> = {};
46105 for (const name of wantVars) {
src/harness/runner.tsmodified+1−5View file
@@ -14,8 +14,6 @@ export interface MatlabSources {
1414 buildProblem: string;
1515 /** laplace2d_bdata.m source */
1616 bdata: string;
17- /** laplace2d_bdata_branch.m source */
18- bdataBranch: string;
1917 /** solver.m source of the solver under test */
2018 solver: string;
2119 }
@@ -50,10 +48,9 @@ export function runPoint(req: RunPointRequest): RunPoint {
5048 const { instance, n } = req;
5149 const repeats = req.repeats ?? 3;
5250 const wantGrid = req.wantGrid ?? false;
53- const family = instance.family === "branch-point" ? 1 : 0;
5451 const main = [
5552 "% generated by the fastandaccurate harness",
56- `prob = build_problem(${numLiteral(instance.a)}, ${numLiteral(instance.k)}, ${numLiteral(instance.d)}, ${family}, ${wantGrid ? 1 : 0});`,
53+ `prob = build_problem(${numLiteral(instance.a)}, ${numLiteral(instance.k)}, ${numLiteral(instance.d)}, ${wantGrid ? 1 : 0});`,
5754 `tic; out = solver(prob, ${n}); res_cold = toc;`,
5855 `res_times = zeros(${repeats}, 1);`,
5956 `for irep = 1:${repeats}`,
@@ -68,7 +65,6 @@ export function runPoint(req: RunPointRequest): RunPoint {
6865 {
6966 "build_problem.m": req.sources.buildProblem,
7067 "laplace2d_bdata.m": req.sources.bdata,
71- "laplace2d_bdata_branch.m": req.sources.bdataBranch,
7268 "solver.m": req.sources.solver,
7369 },
7470 ["res_cold", "res_times", "res_ueval", "res_ugrid"]
src/problems/laplace2d/exact.tsmodified+12−42View file
@@ -22,56 +22,26 @@ export function boundaryPoint(inst: Laplace2dInstance, t: number) {
2222 return { x: r * Math.cos(t), y: r * Math.sin(t) };
2323 }
2424
25-/** Boundary point at parameter phi pushed a distance d along the outward
26- * unit normal. */
27-function pushedPoint(inst: Laplace2dInstance, phi: number, d: number) {
28- const { a, k } = inst;
29- const r = 1 + a * Math.cos(k * phi);
30- const dr = -a * k * Math.sin(k * phi);
31- const bx = r * Math.cos(phi);
32- const by = r * Math.sin(phi);
33- const dx = dr * Math.cos(phi) - r * Math.sin(phi);
34- const dy = dr * Math.sin(phi) + r * Math.cos(phi);
35- const sp = Math.hypot(dx, dy);
36- return { x: bx + (d * dy) / sp, y: by - (d * dx) / sp };
37-}
38-
39-/** The three log-family sources: boundary points at phi_j pushed a
25+/** The three exact-solution sources: boundary points at phi_j pushed a
4026 * distance d along the outward normal, strengths [1.0, -0.6, 0.8]. */
4127 export function sources(inst: Laplace2dInstance): Source[] {
28+ const { a, k, d } = inst;
4229 const strengths = [1.0, -0.6, 0.8];
4330 return strengths.map((c, j) => {
44- const p = pushedPoint(inst, (2 * Math.PI * j) / 3 + 0.4, inst.d);
45- return { ...p, c };
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 };
4640 });
4741 }
4842
49-/** The branch-point family's singularity: the boundary point at
50- * parameter 0.4 pushed d along the outward normal, with its polar
51- * angle (the branch cut runs radially outward from it). */
52-export function branchPoint(inst: Laplace2dInstance) {
53- const p = pushedPoint(inst, 0.4, inst.d);
54- return { ...p, theta0: Math.atan2(p.y, p.x) };
55-}
56-
57-/** The exact solution's singular points, for display. */
58-export function singularities(inst: Laplace2dInstance): { x: number; y: number }[] {
59- return inst.family === "branch-point" ? [branchPoint(inst)] : sources(inst);
60-}
61-
62-/** Exact solution. Log family: u = sum_j c_j log|x - s_j|. Branch
63- * family: u = Re sqrt(w), w = -(z - z0) e^{-i theta0}, evaluated in
64- * real arithmetic as sqrt((|w| + Re w)/2); the branch cut is the radial
65- * ray from z0 away from the origin, outside the domain. */
43+/** Exact solution u(x, y) = sum_j c_j log|x - s_j|. */
6644 export function exactU(inst: Laplace2dInstance, x: number, y: number): number {
67- if (inst.family === "branch-point") {
68- const b = branchPoint(inst);
69- const dx = x - b.x;
70- const dy = y - b.y;
71- const wre = -(dx * Math.cos(b.theta0) + dy * Math.sin(b.theta0));
72- const wim = -(dy * Math.cos(b.theta0) - dx * Math.sin(b.theta0));
73- return Math.sqrt((Math.hypot(wre, wim) + wre) / 2);
74- }
7545 let u = 0;
7646 for (const s of sources(inst)) {
7747 u += s.c * 0.5 * Math.log((x - s.x) ** 2 + (y - s.y) ** 2);
src/problems/laplace2d/matlab/build_problem.mmodified+9−24View file
@@ -1,23 +1,17 @@
1-function prob = build_problem(a, k, d, family, wantGrid)
1+function prob = build_problem(a, k, d, wantGrid)
22 % BUILD_PROBLEM Assemble the problem struct for laplace-dirichlet-2d.
33 %
44 % The domain is the star-shaped region bounded by
55 % x(t) = r(t) [cos t; sin t], r(t) = 1 + a cos(k t), t in [0, 2 pi).
6-% The exact solution depends on the data family:
7-% family 0 (log-sources): u(x) = sum_j c_j log|x - s_j|, with three
8-% point sources s_j outside the domain: s_j is the boundary point at
9-% parameter phi_j = 2 pi (j-1)/3 + 0.4 pushed a distance d along the
10-% outward unit normal, with strengths c = [1.0; -0.6; 0.8].
11-% family 1 (branch-point): u(z) = Re sqrt(-(z - z0) e^{-i th0}), with
12-% z0 the boundary point at parameter 0.4 pushed a distance d along
13-% the outward unit normal and th0 its polar angle (see
14-% laplace2d_bdata_branch.m).
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].
1510 %
1611 % The solver receives only the curve (with derivatives), the Dirichlet
1712 % data g as a function of the boundary parameter, and the points where
18-% the solution is requested. The singularities exist here only to
19-% manufacture the data; a submitted solver must not use knowledge of
20-% them.
13+% the solution is requested. The sources exist here only to manufacture
14+% the data; a submitted solver must not use knowledge of them.
2115 %
2216 % Fields of prob:
2317 % curve @(t) -> [x y] boundary point, t column vector
@@ -28,11 +22,7 @@ function prob = build_problem(a, k, d, family, wantGrid)
2822 % vizXY m x 2 grid points where uGrid is requested
2923 % (m = 0 when no visualization is wanted)
3024
31-if family == 1
32- phi = 0.4;
33-else
34- phi = 2*pi*[0; 1; 2]/3 + 0.4;
35-end
25+phi = 2*pi*[0; 1; 2]/3 + 0.4;
3626 c = [1.0; -0.6; 0.8];
3727 rphi = 1 + a*cos(k*phi);
3828 bx = rphi.*cos(phi);
@@ -49,12 +39,7 @@ prob.curveD = @(t) [-a*k*sin(k*t).*cos(t) - (1 + a*cos(k*t)).*sin(t), ...
4939 -a*k*sin(k*t).*sin(t) + (1 + a*cos(k*t)).*cos(t)];
5040 prob.curveDD = @(t) [(-a*k*k*cos(k*t) - 1 - a*cos(k*t)).*cos(t) + 2*a*k*sin(k*t).*sin(t), ...
5141 (-a*k*k*cos(k*t) - 1 - a*cos(k*t)).*sin(t) - 2*a*k*sin(k*t).*cos(t)];
52-if family == 1
53- th0 = atan2(sy(1), sx(1));
54- prob.g = @(t) laplace2d_bdata_branch(t, a, k, sx(1), sy(1), th0);
55-else
56- prob.g = @(t) laplace2d_bdata(t, a, k, sx, sy, c);
57-end
42+prob.g = @(t) laplace2d_bdata(t, a, k, sx, sy, c);
5843
5944 % Evaluation points: 32 rays, radial fractions 0.1..0.9, plus the origin
6045 % (289 points). The rule must match evalPoints() in
src/problems/laplace2d/matlab/laplace2d_bdata_branch.mdeleted+0−20View file
@@ -1,20 +0,0 @@
1-function g = laplace2d_bdata_branch(t, a, k, x0, y0, th0)
2-% LAPLACE2D_BDATA_BRANCH Dirichlet data for the branch-point family.
3-% The exact solution is u*(z) = Re sqrt(w), w = -(z - z0) e^{-i th0},
4-% with z0 = x0 + i y0 the branch point and th0 its polar angle. The
5-% principal branch cut w <= 0 maps to the radial ray from z0 away from
6-% the origin, which never meets a domain star-shaped about the origin,
7-% so u* is harmonic on the closed domain. Evaluated in real arithmetic:
8-% Re sqrt(w) = sqrt((|w| + Re w)/2).
9-% Called through the prob.g handle built in build_problem; solvers see
10-% only that handle. Must match the branch-point case of exactU in
11-% src/problems/laplace2d/exact.ts.
12-r = 1 + a*cos(k*t);
13-x = r.*cos(t);
14-y = r.*sin(t);
15-dx = x - x0;
16-dy = y - y0;
17-wre = -(dx*cos(th0) + dy*sin(th0));
18-wim = -(dy*cos(th0) - dx*sin(th0));
19-g = sqrt((sqrt(wre.^2 + wim.^2) + wre)/2);
20-end
src/problems/laplace2d/spec.tsmodified+1−23View file
@@ -7,9 +7,6 @@
77 export const PROBLEM_ID = "laplace-dirichlet-2d";
88 export const PROBLEM_VERSION = 1;
99
10-/** The family the exact solution is manufactured from. */
11-export type DataFamily = "log-sources" | "branch-point";
12-
1310 export interface Laplace2dInstance {
1411 /** Short stable identifier used in results and URLs. */
1512 id: string;
@@ -17,9 +14,8 @@ export interface Laplace2dInstance {
1714 /** Boundary r(t) = 1 + a cos(k t). */
1815 a: number;
1916 k: number;
20- /** Distance of the exact solution's singularities beyond the boundary. */
17+ /** Distance of the exact solution's sources beyond the boundary. */
2118 d: number;
22- family: DataFamily;
2319 description: string;
2420 }
2521
@@ -30,7 +26,6 @@ export const INSTANCES: Laplace2dInstance[] = [
3026 a: 0,
3127 k: 0,
3228 d: 0.5,
33- family: "log-sources",
3429 description:
3530 "The unit disk with sources half a radius beyond the boundary. " +
3631 "Every reasonable method should reach high accuracy quickly.",
@@ -41,7 +36,6 @@ export const INSTANCES: Laplace2dInstance[] = [
4136 a: 0.2,
4237 k: 3,
4338 d: 0.4,
44- family: "log-sources",
4539 description:
4640 "A gently star-shaped domain; the data continues comfortably past " +
4741 "the boundary, so geometric convergence is attainable but the " +
@@ -53,26 +47,11 @@ export const INSTANCES: Laplace2dInstance[] = [
5347 a: 0.3,
5448 k: 5,
5549 d: 0.08,
56- family: "log-sources",
5750 description:
5851 "A wavier domain with sources only 0.08 beyond the boundary. The " +
5952 "data barely continues past the boundary, which defeats methods " +
6053 "whose representation assumes it does.",
6154 },
62- {
63- id: "star-branch",
64- label: "3-lobe star, branch-point data",
65- a: 0.2,
66- k: 3,
67- d: 0.4,
68- family: "branch-point",
69- description:
70- "Identical geometry and singularity distance to star-medium, but " +
71- "the data comes from a branch-point singularity (the real part of " +
72- "a complex square root) rather than log point charges, so the pair " +
73- "tests whether behavior depends on the data family or only on its " +
74- "singularity distance.",
75- },
7655 ];
7756
7857 export function getInstance(id: string): Laplace2dInstance {
@@ -90,7 +69,6 @@ export function canonicalSpec(inst: Laplace2dInstance) {
9069 return {
9170 a: inst.a,
9271 d: inst.d,
93- family: inst.family,
9472 instance: inst.id,
9573 k: inst.k,
9674 problem: PROBLEM_ID,
src/solvers/chunkie-dlp/solver.madded+64−0View file
@@ -0,0 +1,64 @@
1+function out = solver(prob, n)
2+% chunkie double-layer BIE for the interior Dirichlet Laplace problem.
3+%
4+% The same second-kind double-layer formulation as nystrom-dlp,
5+% (D - I/2) sigma = g, but discretized and solved by
6+% chunkie (https://github.com/fastalgorithms/chunkie), a production
7+% MATLAB toolbox for boundary integral equations in 2D: the curve is
8+% panelized into n uniform 16th-order Gauss-Legendre chunks, chunkermat
9+% assembles the system with high-order singular quadrature, the dense
10+% system is solved directly, and chunkerkerneval evaluates the potential
11+% with corrected quadrature for targets near the boundary. The package
12+% is fetched by mip on first use.
13+%
14+% n : number of chunks (16 points each).
15+
16+mip load --install magland/magland/chunkie;
17+
18+chnkr = chunkerfuncuni(@(t) fcurve(t, prob), n);
19+
20+% Dirichlet data at the nodes: for this problem family the curve
21+% parameter is the polar angle, so it is recovered from the node
22+% coordinates.
23+xy = chnkr.r(:, :);
24+tt = mod(atan2(xy(2, :), xy(1, :)), 2*pi);
25+rhs = prob.g(tt(:));
26+
27+fkern = kernel('lap', 'd');
28+sysmat = chunkermat(chnkr, fkern);
29+sysmat = sysmat - 0.5*eye(chnkr.npt);
30+sigma = sysmat \ rhs;
31+
32+out = struct();
33+out.uEval = eval_targets(chnkr, fkern, sigma, prob.evalXY);
34+if size(prob.vizXY, 1) > 0
35+ out.uGrid = eval_targets(chnkr, fkern, sigma, prob.vizXY);
36+else
37+ out.uGrid = zeros(0, 1);
38+end
39+
40+end
41+
42+function u = eval_targets(chnkr, fkern, sigma, XY)
43+% Direct (unaccelerated) evaluation, in blocks to bound memory. accel is
44+% disabled because chunkie's FMM acceleration binds to the fmm2d
45+% library, which is not available in this embedded numbl runtime; at
46+% these sizes direct evaluation is cheap anyway.
47+opts = struct();
48+opts.accel = false;
49+m = size(XY, 1);
50+u = zeros(m, 1);
51+B = 2000;
52+for i0 = 1:B:m
53+ i1 = min(i0 + B - 1, m);
54+ ub = chunkerkerneval(chnkr, fkern, sigma, XY(i0:i1, :).', opts);
55+ u(i0:i1) = ub(:);
56+end
57+end
58+
59+function [r, d, d2] = fcurve(t, prob)
60+tt = t(:);
61+r = prob.curve(tt).';
62+d = prob.curveD(tt).';
63+d2 = prob.curveDD(tt).';
64+end
src/solvers/index.tsmodified+18−0View file
@@ -45,6 +45,24 @@ export const SOLVERS: SolverManifest[] = [
4545 backend: "cpu",
4646 sweepN: [16, 24, 32, 48, 64, 96, 128, 192, 256, 384, 512, 768],
4747 },
48+ {
49+ id: "chunkie-dlp",
50+ name: "chunkie double-layer BIE",
51+ description:
52+ "The same second-kind double-layer formulation as nystrom-dlp, " +
53+ "discretized by chunkie, a production MATLAB boundary-integral " +
54+ "toolbox: n uniform 16th-order Gauss-Legendre panels, high-order " +
55+ "singular quadrature in the assembly, a direct dense solve, and " +
56+ "near-corrected evaluation of the potential. chunkie's default " +
57+ "quadrature tolerances cap the attainable accuracy near 1e-11, and " +
58+ "the timings include the cost of running a general-purpose library " +
59+ "through numbl. The package is fetched by mip on first use, so the " +
60+ "first run in a session spends tens of seconds downloading it; " +
61+ "later runs do not.",
62+ version: "1.0.0",
63+ backend: "cpu",
64+ sweepN: [2, 3, 4, 6, 8, 12, 16, 24, 32, 48],
65+ },
4866 ];
4967
5068 export function getSolver(id: string): SolverManifest {
test/solver-test.tsmodified+10−13View file
@@ -9,6 +9,12 @@ import { INSTANCES, getInstance } from "../src/problems/laplace2d/spec";
99 import { SOLVERS } from "../src/solvers";
1010 import { runPoint, type MatlabSources } from "../src/harness/runner";
1111 import { runSweep } from "../src/harness/sweep";
12+import { setNumblFileIO } from "../src/harness/numblRun";
13+import { NodeFileIOAdapter } from "../src/cli/nodeFileIO";
14+
15+// Solvers that mip-install packages (chunkie-dlp) need file I/O; in node
16+// that is the curl-backed adapter.
17+setNumblFileIO((vfs) => new NodeFileIOAdapter(vfs));
1218
1319 const root = join(dirname(fileURLToPath(import.meta.url)), "..");
1420 const read = (p: string) => readFileSync(join(root, p), "utf-8");
@@ -16,11 +22,11 @@ const read = (p: string) => readFileSync(join(root, p), "utf-8");
1622 const base = {
1723 buildProblem: read("src/problems/laplace2d/matlab/build_problem.m"),
1824 bdata: read("src/problems/laplace2d/matlab/laplace2d_bdata.m"),
19- bdataBranch: read("src/problems/laplace2d/matlab/laplace2d_bdata_branch.m"),
2025 };
2126 const solverSources: Record<string, MatlabSources> = {
2227 mfs: { ...base, solver: read("src/solvers/mfs/solver.m") },
2328 "nystrom-dlp": { ...base, solver: read("src/solvers/nystrom-dlp/solver.m") },
29+ "chunkie-dlp": { ...base, solver: read("src/solvers/chunkie-dlp/solver.m") },
2430 };
2531
2632 // Best relMax each solver must reach over its full sweep. On star-hard,
@@ -28,18 +34,9 @@ const solverSources: Record<string, MatlabSources> = {
2834 // beyond the data's singularities there, and if it suddenly reached high
2935 // accuracy the instance would no longer be testing what the spec says.
3036 const mustReach: Record<string, Record<string, number>> = {
31- mfs: {
32- "disk-easy": 1e-12,
33- "star-medium": 1e-12,
34- "star-hard": 1e-2,
35- "star-branch": 1e-10,
36- },
37- "nystrom-dlp": {
38- "disk-easy": 1e-10,
39- "star-medium": 1e-10,
40- "star-hard": 1e-8,
41- "star-branch": 1e-10,
42- },
37+ mfs: { "disk-easy": 1e-12, "star-medium": 1e-12, "star-hard": 1e-2 },
38+ "nystrom-dlp": { "disk-easy": 1e-10, "star-medium": 1e-10, "star-hard": 1e-8 },
39+ "chunkie-dlp": { "disk-easy": 1e-10, "star-medium": 1e-10, "star-hard": 1e-9 },
4340 };
4441 const mustNotReach: Record<string, Record<string, number>> = {
4542 mfs: { "star-hard": 1e-8 },