Three new instances, near-field target sets, and an adaptive timing policy
flower-15 keeps the difficulty in the geometry alone: fifteen lobes, a
curvature radius of 0.014, sources a comfortable 0.25 out. square-corners
moves it from everywhere to somewhere, on a second boundary family (the
superellipse |x|^p + |y|^p = 1, analytic for even p, four corners of
radius about 1.4/p). star-nearfield changes the question rather than the
problem: star-hard's domain and data, with 32 targets 0.005 to 0.05 inside
the boundary along the inward normal, so that close evaluation is measured
by itself. There chunkie-dlp reaches 1.8e-11, where nystrom-dlp's
uncorrected quadrature is held at 1.2e-2 and mfs at its conditioning wall
near 2e-6.
The evaluation-point count is therefore per instance (289, 305 or 321) and
a solver must read it from prob.evalXY. Sweeps get per-instance resolution
lists, since a harder instance needs a different range and putting those
resolutions in every sweep would only slow the easy ones.
Timing becomes adaptive (src/harness/timing.ts, shared by the numbl and
MATLAB runners): timed runs continue until they have used a time budget,
subject to a floor and a cap, which is what makes the cheap end of the
curves reproducible. The result protocol records the policy.
chunkie now comes from mip, whose fmm2d package ships a compiled MEX
binary, so chunkie-dlp (3.0.0) runs with chunkie's own defaults and its
FMM acceleration rather than with acceleration deliberately disabled.
Solver expectations move to test/expected.ts, keyed by source directory,
so that a -mat entry and its numbl twin are held to the same accuracy.
25 changed files+1203−278
README.mdmodified+19−6View file
@@ -17,9 +17,11 @@ Solvers are MATLAB function files. Most run via
1717 [numbl](https://numbl.org) (MATLAB syntax in the browser and in node),
1818 both on the site and from the command line; some run only in real
1919 MATLAB through the command line, and their results are marked as not
20-reproducible in the browser. Each problem defines its own interface and
21-instances in a written specification; interfaces are per problem rather
22-than shared.
20+reproducible in the browser. Two registry entries may share one file:
21+the `-mat` solvers are their numbl twin's `solver.m` run in real MATLAB,
22+so that pair of curves measures the runtime rather than the method. Each
23+problem defines its own interface and instances in a written
24+specification; interfaces are per problem rather than shared.
2325
2426 ## Problems
2527
@@ -49,9 +51,20 @@ npx https://concept-collection.github.io/fastandaccurate/cli.tgz run --label "my
4951 Note that npx caches by the exact URL string; the site offers the URL
5052 with a `?v=<commit>` suffix so each deployment is a fresh install.
5153
52-Useful flags: `--instance <id>`, `--solver <id>`, `--repeats N`,
53-`--max-n N`, `--out dir`. To benchmark your own solver, point the
54-harness at a MATLAB function file implementing the problem's interface:
54+The solvers whose runtime is `matlab` need `matlab` on the PATH; the run
55+skips them when it is absent. `chunkie-dlp` needs one thing more, the
56+[mip](https://mip.sh) package manager on the MATLAB path, from which the
57+harness installs chunkie and its FLAM and fmm2d dependencies on first
58+use. Taking chunkie from mip rather than from a source clone is what
59+makes its accelerated code path available without a Fortran compiler on
60+the machine, since the mip fmm2d package ships a compiled MEX binary per
61+platform.
62+
63+Useful flags: `--instance <id>`, `--solver <id>`, `--repeats N` (the
64+minimum timed runs per point; each point is then repeated until it has
65+used the `--time-budget`, 0.5 s by default), `--max-n N`, `--out dir`.
66+To benchmark your own solver, point the harness at a MATLAB function file
67+implementing the problem's interface:
5568
5669 ```
5770 npx https://concept-collection.github.io/fastandaccurate/cli.tgz run \
docs/problems/laplace-dirichlet-2d.mdmodified+157−27View file
@@ -10,11 +10,24 @@ Solve the Laplace equation
1010
1111 Δu = 0 in Ω, u = g on ∂Ω,
1212
13-where Ω is the star-shaped plane domain bounded by the curve
13+where Ω is the plane domain bounded by the curve
1414
15- x(t) = r(t) (cos t, sin t), r(t) = 1 + a cos(k t), t ∈ [0, 2π),
15+ x(t) = r(t) (cos t, sin t), t ∈ [0, 2π),
1616
17-and g is Dirichlet data specified below. The problem is posed in the
17+and g is Dirichlet data specified below. Two families of radius function
18+are used, both analytic and both parameterized by the polar angle:
19+
20+ star: r(t) = 1 + a cos(k t)
21+ rounded square: r(t) = (cos^p t + sin^p t)^(−1/p), p even
22+
23+The second traces the superellipse |x|^p + |y|^p = 1: four sides, straight
24+to within about 1/p², meeting at four corners on the diagonals whose
25+radius of curvature is about 1.4/p. Since p is even, cos^p t is a
26+polynomial in cos t and the curve is analytic, so `curveDD` exists
27+everywhere and a solver sees the same interface as for a star. A genuine
28+corner, where the curve is only piecewise smooth and the layer density is
29+singular, is a different problem and not this one; what this family gives
30+is an arbitrarily *near* corner. The problem is posed in the
1831 continuum: a solver receives the curve and the data as functions, chooses
1932 its own discretization, and is scored on the values it returns at a fixed
2033 set of evaluation points.
@@ -38,11 +51,58 @@ not use knowledge of the sources; they exist only to manufacture g.
3851
3952 ## Official instances
4053
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 |
54+| id | shape | a | k | p | d | character |
55+|---|---|---|---|---|---|---|
56+| disk-easy | star | 0 | 0 | | 0.5 | the unit disk, distant sources |
57+| star-medium | star | 0.2 | 3 | | 0.4 | mild geometry, comfortable continuation |
58+| star-hard | star | 0.3 | 5 | | 0.08 | wavy boundary, data barely continues |
59+| flower-15 | star | 0.2 | 15 | | 0.25 | fine features everywhere, data continues |
60+| square-corners | rounded square | | | 100 | 0.25 | near corners, targets inside them |
61+| star-nearfield | star | 0.3 | 5 | | 0.08 | star-hard, with targets against the boundary |
62+
63+The hard instances are hard in different ways, which is the point of
64+having several. On star-hard the geometry is mild and the data barely
65+continues past the boundary; on flower-15 the data continues comfortably
66+(the sources really are 0.25 from the boundary, not merely 0.25 along a
67+normal) and the boundary itself is the difficulty: fifteen lobes, with
68+curvature up to 69, a curvature radius of 0.014 against a domain of size
69+1. Resolving it costs a few hundred boundary nodes before any method
70+converges at all, and the fixed evaluation rule below puts the closest
71+evaluation point 0.0295 from the boundary, so a method whose quadrature
72+degrades near the boundary is held back further. Note that k = 15 is
73+divisible by 3, so the three sources sit at equivalent phases of the lobe
74+pattern; this is deliberate, and keeps the instance from favoring any
75+particular angular sector.
76+
77+square-corners moves the difficulty from everywhere to somewhere. Its
78+corners have radius 0.0142, one percent of the side, and its sides are
79+straight to a thousandth: four short stretches of boundary demand all the
80+resolution and the rest demands almost none. A method that can grade its
81+discretization has something to exploit here that the star instances never
82+offer, and the extra evaluation points described below sit inside the
83+corners, at 0.005 to 0.05 from the boundary, so that near-field accuracy
84+is measured where the geometry is worst rather than only where it is
85+convenient.
86+
87+star-nearfield changes the question rather than the problem. Its domain
88+and its data are exactly star-hard's, so the equation a solver has to
89+solve is the same one; what differs is where the answer is wanted. Besides
90+the standard points it carries 32 targets on the inward normal at eight
91+parameters around the curve, at the same four distances the near-corner
92+set uses, so the closest sits 0.005 from a boundary of perimeter 9.02. A
93+method that evaluates its own representation with an uncorrected
94+quadrature rule loses accuracy there for reasons that have nothing to do
95+with how well it solved the equation: for the periodic trapezoid rule the
96+error at a target a distance δ inside behaves like exp(−2πδ/h) in the node
97+spacing h, so at n = 1024 nodes (h = 0.0088) one and the same computed
98+density is accurate to 6e−14 at δ = 0.05, 6.6e−6 at 0.02, 2.7e−3 at 0.01
99+and only 6.3e−2 at 0.005. The instance therefore measures near-field
100+evaluation by itself, and the pair (star-hard, star-nearfield) isolates
101+it: identical geometry, identical data, identical discretizations, one
102+question changed. Note that a method whose representation is smooth up to
103+the boundary, such as the MFS with its charges outside the domain, pays
104+nothing here; the difficulty belongs to the representation and not to the
105+problem.
46106
47107 Committed results exist only at these instances, so every solver is
48108 compared on identical inputs. The parameters may be varied freely in the
@@ -50,34 +110,79 @@ site's interactive views, but such runs are not recorded.
50110
51111 ## What is scored
52112
53-Each instance fixes 289 evaluation points: 32 rays at angles
113+Each instance fixes at least 289 evaluation points: 32 rays at angles
54114 θ_j = 2πj/32 + 0.13 (j = 0, …, 31), radial fractions
55115 ρ ∈ {0.1, 0.2, …, 0.9} along each ray (the point at fraction ρ on ray
56116 θ is ρ r(θ) (cos θ, sin θ)), ordered radius-major, plus the origin last.
117+
118+An instance with corners adds a **near-corner set** after those: for each
119+corner in order (the diagonals θ_c = π/4 + jπ/2, j = 0, …, 3), the four
120+points at distances δ ∈ {0.005, 0.01, 0.02, 0.05} inside the boundary
121+along that diagonal, at radius r(θ_c) − δ. On square-corners this gives
122+16 further points and 305 in all, and their true perpendicular distances
123+to the boundary are 0.005, 0.010, 0.019 and 0.042, the smallest of them a
124+third of the corner's own radius of curvature.
125+
126+An instance carrying the **near-boundary set** adds, after any near-corner
127+set, four points per parameter: for each of the eight parameters
128+t_i = 2πi/8 + 0.07 (i = 0, …, 7), the four points at the same distances
129+δ ∈ {0.005, 0.01, 0.02, 0.05} inside the boundary along the inward *unit
130+normal* at t_i, ordered with the parameter outer and the distance inner.
131+On star-nearfield this gives 32 further points and 321 in all. Because
132+every δ stays below the smallest radius of curvature
133+there (0.0721), the distance from such a point to the curve is exactly δ,
134+unlike the near-corner set whose radial offset makes the true distances a
135+little smaller. The eight parameters are offset by 0.07 so that they land
136+at no special phase of the lobe pattern, and 8 is coprime to k = 5, so
137+they sample crests, flanks and valleys alike.
138+
139+The count is therefore per instance, and a solver must return one value
140+per row of `prob.evalXY` rather than assume 289.
141+
57142 The solver returns u at exactly these points. Reported errors are
58143
59144 relMax = max_i |u_i − u*_i| / max_i |u*_i|,
60145 relL2 = ( Σ_i (u_i − u*_i)² / Σ_i u*_i² )^{1/2}.
61146
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.
147+relMax is the headline number; both are recorded. Note that the fraction
148+0.9 bounds the *radial* distance, not the distance to the boundary. On a
149+strongly wiggly instance the boundary runs nearly radially along the
150+flanks of the lobes, so the perpendicular distance from an evaluation
151+point to it is several times smaller than the radial one: 0.0295 at the
152+closest point on flower-15, against 0.0464 on star-hard and 0.077 on
153+star-medium. Close evaluation therefore enters on the wiggly instances,
154+though not in its severe form. The severe form is star-nearfield, whose
155+near-boundary set puts targets a fraction of a node spacing from the
156+boundary. Adding it needed no version bump: it is a new instance, and no
157+number measured on an existing one changes.
69158
70159 ## Timing protocol
71160
72161 One run is a full call of the solver, including its own discretization,
73162 assembly, solve, and evaluation. The harness performs two untimed warmup
74-runs (the first of which is recorded separately as the cold time), then N
75-timed runs (N = 5 unless stated), and reports the **fastest** of them as
76-the solve time. The minimum is used rather than the mean or median
77-because everything that interferes with a measurement (scheduling, other
78-load, residual JIT compilation) only ever adds time, so the fastest run
79-is the least contaminated estimate of the solver's own cost; every
80-individual timing is recorded in the result file regardless. Where a
163+runs (the first of which is recorded separately as the cold time), then
164+timed runs, and reports the **fastest** of them as the solve time. The
165+minimum is used rather than the mean or median because everything that
166+interferes with a measurement (scheduling, other load, residual JIT
167+compilation) only ever adds time, so the fastest run is the least
168+contaminated estimate of the solver's own cost; every individual timing
169+is recorded in the result file regardless.
170+
171+The number of timed runs is not fixed, because a fixed count times a
172+cheap point badly: a solve of 0.1 ms is dominated by whatever else the
173+machine was doing during those few hundred microseconds, and the minimum
174+of five such samples still scatters by a factor of several from one sweep
175+to the next, which made the low-resolution end of the curves noisy. So
176+timed runs continue until they have accumulated a **time budget** (0.5 s
177+by default), subject to a floor on the count (5) and a cap (50 runs).
178+A solve of a fraction of a millisecond is therefore measured 50 times and
179+an expensive one 5 times, and each point is sampled enough for its
180+minimum to mean something. The cap also bounds how many timings a result
181+file carries, since it records all of them. The policy is recorded in the result file as
182+`protocol.minTimedRuns`, `protocol.timeBudgetSeconds`, and
183+`protocol.maxTimedRuns`; how many runs a given point actually took is the
184+length of its `solveSecondsAll`. Results written before this policy
185+record a fixed `protocol.timedRuns` instead. Where a
81186 sweep shares one process across resolutions, as the MATLAB runner does,
82187 warmup runs precede the sweep as well, so that the first resolution does
83188 not absorb the session's one-time costs. All timing is tic/toc inside the
@@ -101,7 +206,8 @@ A solver is a MATLAB function file
101206
102207 where n is the solver's own resolution parameter (its meaning is the
103208 solver's choice; the standard sweep list is declared in the solver's
104-manifest) and prob is a struct with fields
209+manifest, which may also declare a different list for an instance that
210+needs a different range of resolutions) and prob is a struct with fields
105211
106212 | field | meaning |
107213 |---|---|
@@ -109,10 +215,10 @@ manifest) and prob is a struct with fields
109215 | `prob.curveD` | first derivative of the curve with respect to t |
110216 | `prob.curveDD` | second derivative |
111217 | `prob.g` | `@(t) -> g`, Dirichlet data at boundary parameter t |
112-| `prob.evalXY` | 289×2, the evaluation points |
218+| `prob.evalXY` | m×2, the evaluation points (289; 305 with corners, 321 with the near-boundary set) |
113219 | `prob.vizXY` | m×2, visualization grid points (m = 0 when not requested) |
114220
115-The return value is a struct: `out.uEval` (289×1, required) and
221+The return value is a struct: `out.uEval` (m×1, required) and
116222 `out.uGrid` (m×1; `[]` when `prob.vizXY` is empty). The solver must not
117223 reconstruct the sources analytically or otherwise special-case the known
118224 solution; submissions are reviewed for this.
@@ -120,7 +226,7 @@ solution; submissions are reviewed for this.
120226 ## Visualization grid
121227
122228 When requested, `prob.vizXY` lists a 200×200 grid of points over the
123-bounding square [−R, R]², R = 1.05(1 + |a|), with flat index
229+bounding square [−R, R]², R = 1.05 max_t r(t), with flat index
124230 p = ix·200 + iy for x = xs[ix], y = xs[iy] (y varies fastest, MATLAB
125231 meshgrid column order). Points outside Ω are included and the viewer
126232 masks them; grid values are never scored.
@@ -132,3 +238,27 @@ problem does not test corner handling, nonsmooth data, or interior
132238 sources (a nonzero right-hand side would exclude plain boundary-integral
133239 methods; that belongs to a different problem). The domain family is
134240 star-shaped by construction, which some methods can exploit.
241+
242+Every boundary here is analytic, which bounds what the measured curves
243+can show, and the instances differ in how much that bound bites. On the
244+star instances the features are of one size everywhere: flower-15 is
245+fifteen copies of the same lobe, so a method that refines adaptively has
246+nothing to refine toward, and a global spectral rule, which must resolve
247+the finest feature everywhere, loses nothing by doing so. There the
248+geometric convergence of the global trapezoid rule is hard to beat.
249+square-corners is the deliberate counterweight, with one localized feature
250+and targets inside it, and the ranking there is different in kind rather
251+than in degree. star-nearfield is a second counterweight, aimed at
252+evaluation rather than at geometry: what orders the methods there is
253+whether they correct their quadrature near the boundary at all, and a
254+method that does not is held at 1e−2 where one that does reaches 1e−11.
255+
256+What remains untested is a genuine corner. The rounded square is smooth,
257+so its layer density is smooth; a boundary with an actual vertex has a
258+density that behaves like a power of the distance to it, needs graded
259+refinement or a corner-specific quadrature, and can turn a well-posed
260+second-kind equation into something that needs care. That is a different
261+problem, not a harder instance of this one. Also untested is whether
262+adaptive refinement pays: the panel solver here places its panels
263+uniformly, so on square-corners it resolves the corner by refining the
264+flat sides at the same time, and the benefit of grading is not measured.
src/app/components/DomainView.tsxmodified+11−6View file
@@ -6,6 +6,7 @@ import type { Laplace2dInstance } from "../../problems/laplace2d/spec";
66 import {
77 boundaryPoint,
88 evalPoints,
9+ maxRadius,
910 sources,
1011 } from "../../problems/laplace2d/exact";
1112
@@ -15,6 +16,9 @@ function token(name: string): string {
1516
1617 export function DomainView({ inst }: { inst: Laplace2dInstance }) {
1718 const canvasRef = useRef<HTMLCanvasElement>(null);
19+ // The count varies by instance: 289 normally, and more where a
20+ // near-corner or near-boundary set is added.
21+ const nEval = evalPoints(inst).length;
1822
1923 useEffect(() => {
2024 const canvas = canvasRef.current;
@@ -31,15 +35,16 @@ export function DomainView({ inst }: { inst: Laplace2dInstance }) {
3135 ctx.scale(dpr, dpr);
3236 ctx.clearRect(0, 0, size, size);
3337
34- const extent = 1 + Math.abs(inst.a) + inst.d + 0.25;
38+ const extent = maxRadius(inst) + inst.d + 0.25;
3539 const s = size / (2 * extent);
3640 const X = (x: number) => size / 2 + x * s;
3741 const Y = (y: number) => size / 2 - y * s;
3842
3943 // domain fill + boundary
4044 ctx.beginPath();
41- for (let i = 0; i <= 512; i++) {
42- const t = (2 * Math.PI * i) / 512;
45+ const nOutline = 4096;
46+ for (let i = 0; i <= nOutline; i++) {
47+ const t = (2 * Math.PI * i) / nOutline;
4348 const p = boundaryPoint(inst, t);
4449 if (i === 0) ctx.moveTo(X(p.x), Y(p.y));
4550 else ctx.lineTo(X(p.x), Y(p.y));
@@ -83,9 +88,9 @@ export function DomainView({ inst }: { inst: Laplace2dInstance }) {
8388 <figure style={{ margin: 0 }}>
8489 <canvas ref={canvasRef} />
8590 <figcaption className="field-caption" style={{ maxWidth: 360 }}>
86- The domain, the evaluation points where solutions are scored
87- (dots), and the exact solution's sources a distance {inst.d} outside
88- the boundary (crosses).
91+ The domain, the {nEval} evaluation points where solutions are
92+ scored (dots), and the exact solution's sources a distance {inst.d}{" "}
93+ outside the boundary (crosses).
8994 </figcaption>
9095 </figure>
9196 );
src/app/components/SolutionSection.tsxmodified+13−6View file
@@ -10,7 +10,7 @@ import {
1010 vizGrid,
1111 VIZ_NGRID,
1212 } from "../../problems/laplace2d/exact";
13-import { SOLVERS, getSolver } from "../../solvers";
13+import { SOLVERS, getSolver, sweepNFor } from "../../solvers";
1414 import { solutionInBrowser } from "../workerClient";
1515 import type { ResultPoint } from "../../harness/resultSchema";
1616 import { FieldView, fieldAbsMax } from "./FieldView";
@@ -70,12 +70,19 @@ export function SolutionSection({ inst }: { inst: Laplace2dInstance }) {
7070
7171 const marks = useMemo(() => evalPoints(inst), [inst]);
7272
73+ // The resolutions offered depend on the instance as well as the solver,
74+ // since a harder geometry sweeps a different range. If the instance
75+ // changes under a selected n that its list does not contain, fall back
76+ // to a resolution partway up the new list.
77+ const ns = useMemo(() => sweepNFor(solver, inst.id), [solver, inst.id]);
78+ const nSel = ns.includes(n) ? n : ns[Math.floor(ns.length * 0.7)];
79+
7380 async function compute() {
7481 setBusy(true);
7582 setError(null);
7683 try {
77- const { point, uGrid } = await solutionInBrowser(inst.id, solverId, n);
78- setComputed({ solverId, n, uGrid, point });
84+ const { point, uGrid } = await solutionInBrowser(inst.id, solverId, nSel);
85+ setComputed({ solverId, n: nSel, uGrid, point });
7986 } catch (err) {
8087 setError(err instanceof Error ? err.message : String(err));
8188 } finally {
@@ -93,7 +100,7 @@ export function SolutionSection({ inst }: { inst: Laplace2dInstance }) {
93100 onChange={(e) => {
94101 const id = e.target.value;
95102 setSolverId(id);
96- const sw = getSolver(id).sweepN;
103+ const sw = sweepNFor(getSolver(id), inst.id);
97104 setN(sw[Math.floor(sw.length * 0.7)]);
98105 }}
99106 >
@@ -106,8 +113,8 @@ export function SolutionSection({ inst }: { inst: Laplace2dInstance }) {
106113 </label>
107114 <label>
108115 n{" "}
109- <select value={n} onChange={(e) => setN(parseInt(e.target.value, 10))}>
110- {solver.sweepN.map((v) => (
116+ <select value={nSel} onChange={(e) => setN(parseInt(e.target.value, 10))}>
117+ {ns.map((v) => (
111118 <option key={v} value={v}>
112119 {v}
113120 </option>
src/app/matlabSources.tsmodified+7−2View file
@@ -1,5 +1,8 @@
11 // 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).
2+// and by the problem page's source listing; the node CLI reads the same
3+// files from disk.
4+
5+import { getSolver, solverSourceDir } from "../solvers";
36
47 const files = import.meta.glob("../{problems,solvers}/**/*.m", {
58 query: "?raw",
@@ -20,6 +23,8 @@ export function matlabBase() {
2023 };
2124 }
2225
26+/** The solver.m of a registry solver, resolved through its manifest so
27+ * that entries sharing a file (mfs and mfs-mat) read the one source. */
2328 export function solverSource(solverId: string): string {
24- return get(`../solvers/${solverId}/solver.m`);
29+ return get(`../solvers/${solverSourceDir(getSolver(solverId))}/solver.m`);
2530 }
src/app/pages/AboutPage.tsxmodified+20−6View file
@@ -31,11 +31,16 @@ export function AboutPage() {
3131 error against compute time, traced out as the solver's resolution
3232 parameter varies. Errors are measured at a fixed set of evaluation
3333 points defined per instance, relative to the reference solution.
34- Timing is two untimed warmup runs followed by several timed runs, of
35- which the fastest is reported: interference only ever adds time, so
36- the fastest run is the least contaminated estimate of the solver's
37- own cost. A run includes the solver's own discretization, assembly,
38- solve, and evaluation. Curves are traced in order of the solver's
34+ Timing is two untimed warmup runs followed by timed runs, of which
35+ the fastest is reported: interference only ever adds time, so the
36+ fastest run is the least contaminated estimate of the solver's own
37+ cost. How many timed runs is not fixed but set by a time budget: a
38+ point keeps repeating until it has used about half a second, with a
39+ floor of five runs and a cap of fifty, so a solve that takes a
40+ fraction of a millisecond is sampled fifty times rather than five,
41+ and its reported time stops depending on what else the machine
42+ happened to be doing. A run includes the solver's own
43+ discretization, assembly, solve, and evaluation. Curves are traced in order of the solver's
3944 resolution parameter, and may double back in time, since a solver's
4045 cost need not increase with resolution.
4146 </p>
@@ -63,12 +68,21 @@ export function AboutPage() {
6368 <p>
6469 This runs the standard sweeps and writes one result JSON per instance
6570 and solver. Useful flags: <code>--instance <id></code>,{" "}
66- <code>--solver <id></code>, <code>--repeats N</code>,{" "}
71+ <code>--solver <id></code>, <code>--repeats N</code> (the
72+ minimum timed runs per point), <code>--time-budget S</code>,{" "}
6773 <code>--max-n N</code>, <code>--out dir</code>. To benchmark your own
6874 solver, point the harness at a MATLAB function file implementing the
6975 problem's solver interface:
7076 </p>
7177 <pre>{`npx ${cliUrl} run --solver-file my_method.m --solver-id my-method`}</pre>
78+ <p>
79+ The solvers that run only in real MATLAB need <code>matlab</code> on
80+ the PATH; the run skips them when it is absent.{" "}
81+ <code>chunkie-dlp</code> needs one thing more, the{" "}
82+ <a href="https://mip.sh">mip</a> package manager on the MATLAB path,
83+ from which the harness installs chunkie and its FLAM and fmm2d
84+ dependencies on first use.
85+ </p>
7286 <p className="small muted">
7387 Note that npx caches by the exact URL string; the <code>?v=</code>{" "}
7488 suffix above ties the command to the current deployment so a later
src/app/pages/ProblemPage.tsxmodified+50−20View file
@@ -5,7 +5,7 @@ import {
55 getInstance,
66 PROBLEM_ID,
77 } from "../../problems/laplace2d/spec";
8-import { SOLVERS } from "../../solvers";
8+import { SOLVERS, solverSourceDir } from "../../solvers";
99 import type { ResultFile, ResultPoint } from "../../harness/resultSchema";
1010 import {
1111 environmentLabel,
@@ -14,6 +14,7 @@ import {
1414 } from "../results";
1515 import { solverColorVar } from "../colors";
1616 import { sweepInBrowser } from "../workerClient";
17+import { DEFAULT_TIMING } from "../../harness/timing";
1718 import { solverSource } from "../matlabSources";
1819 import {
1920 WorkPrecisionChart,
@@ -30,7 +31,7 @@ interface LocalRun {
3031 key: string;
3132 solverId: string;
3233 instanceId: string;
33- repeats: number;
34+ minTimedRuns: number;
3435 points: ResultPoint[];
3536 done: boolean;
3637 }
@@ -44,7 +45,7 @@ export function ProblemPage({ problemId }: { problemId: string }) {
4445 const [hidden, setHidden] = useState<Set<string>>(new Set());
4546 const [running, setRunning] = useState<string | null>(null);
4647 const [runStatus, setRunStatus] = useState<string | null>(null);
47- const [repeats, setRepeats] = useState(3);
48+ const [minRuns, setMinRuns] = useState(DEFAULT_TIMING.minTimedRuns);
4849 const [copied, setCopied] = useState(false);
4950
5051 const inst = getInstance(instanceId);
@@ -118,14 +119,14 @@ export function ProblemPage({ problemId }: { problemId: string }) {
118119 ...rs.filter(
119120 (r) => !(r.solverId === solverId && r.instanceId === instanceId)
120121 ),
121- { key, solverId, instanceId, repeats, points: [], done: false },
122+ { key, solverId, instanceId, minTimedRuns: minRuns, points: [], done: false },
122123 ]);
123124 setRunning(solverId);
124125 try {
125126 const points = await sweepInBrowser(
126127 instanceId,
127128 solverId,
128- repeats,
129+ { ...DEFAULT_TIMING, minTimedRuns: minRuns },
129130 (point, index, total) => {
130131 setRunStatus(
131132 `${solverId} on ${instanceId}: point ${index + 1}/${total} (n = ${point.n}) — rel max error ${point.relMax.toExponential(2)}`
@@ -195,8 +196,11 @@ export function ProblemPage({ problemId }: { problemId: string }) {
195196 Interior Dirichlet Laplace problem on a star-shaped 2D domain
196197 </p>
197198 <p>
198- Solve Δu = 0 on the domain with boundary r(θ) = 1 + a·cos(kθ), with
199- Dirichlet data u = g on the boundary. The data comes from an exact
199+ Solve Δu = 0 on the domain, with Dirichlet data u = g on the
200+ boundary. The boundary is r(θ) = 1 + a·cos(kθ) on most instances,
201+ and on <code>square-corners</code> the superellipse
202+ |x|^p + |y|^p = 1, a square whose corners are rounded to a radius
203+ of about 1.4/p. The data comes from an exact
200204 harmonic function, a sum of three logarithmic point sources placed a
201205 distance d outside the boundary, so errors are measured against the
202206 true solution rather than a reference computation. The distance d
@@ -205,7 +209,11 @@ export function ProblemPage({ problemId }: { problemId: string }) {
205209 representations assume that continuation lose it. A solver receives
206210 the curve (with derivatives), the boundary data as a function of the
207211 boundary parameter, and the evaluation points, and returns solution
208- values at those points. The precise statement, solver interface, and
212+ values at those points. Those points are the same 289 on most
213+ instances; <code>square-corners</code> adds sixteen inside its
214+ corners and <code>star-nearfield</code> thirty-two along the inward
215+ normal, as close as 0.005 to the boundary, where a quadrature rule
216+ with no near-field correction stops converging. The precise statement, solver interface, and
209217 timing protocol are in the <a href={SPEC_URL}>specification</a>.
210218 </p>
211219 <div className="row" style={{ marginTop: 14 }}>
@@ -231,10 +239,21 @@ export function ProblemPage({ problemId }: { problemId: string }) {
231239 <table className="data">
232240 <tbody>
233241 <tr>
234- <th className="left">a</th>
235- <td>{inst.a}</td>
236- <th className="left">k</th>
237- <td>{inst.k}</td>
242+ {inst.shape === "rounded-square" ? (
243+ <>
244+ <th className="left">p</th>
245+ <td>{inst.p}</td>
246+ <th className="left">corners</th>
247+ <td>4</td>
248+ </>
249+ ) : (
250+ <>
251+ <th className="left">a</th>
252+ <td>{inst.a}</td>
253+ <th className="left">k</th>
254+ <td>{inst.k}</td>
255+ </>
256+ )}
238257 <th className="left">d</th>
239258 <td>{inst.d}</td>
240259 </tr>
@@ -247,8 +266,11 @@ export function ProblemPage({ problemId }: { problemId: string }) {
247266 <h2>Solvers</h2>
248267 <p className="small muted" style={{ maxWidth: 640 }}>
249268 Each solver is a MATLAB function file implementing the interface in
250- the <a href={SPEC_URL}>specification</a>; the same file runs in the
251- browser via numbl and from the command line.
269+ the <a href={SPEC_URL}>specification</a>. Most run through numbl, in
270+ the browser and from the command line alike; the rest run only in
271+ real MATLAB, through the command line. The <code>-mat</code> entries
272+ are their numbl twin's file run that way, so each such pair measures
273+ the runtime rather than the method.
252274 </p>
253275 {SOLVERS.map((s) => (
254276 <div
@@ -267,6 +289,7 @@ export function ProblemPage({ problemId }: { problemId: string }) {
267289 {s.runtime === "matlab"
268290 ? "runs in MATLAB via the command line"
269291 : "runs via numbl in the browser and command line"}
292+ {s.sourceDir ? ` · same solver.m as ${s.sourceDir}` : ""}
270293 </span>
271294 </div>
272295 <p className="small" style={{ color: "var(--text-2)" }}>
@@ -282,7 +305,7 @@ export function ProblemPage({ problemId }: { problemId: string }) {
282305 </details>
283306 <a
284307 className="small"
285- href={`${REPO_URL}/blob/main/src/solvers/${s.id}/solver.m`}
308+ href={`${REPO_URL}/blob/main/src/solvers/${solverSourceDir(s)}/solver.m`}
286309 >
287310 view on GitHub
288311 </a>
@@ -332,13 +355,20 @@ export function ProblemPage({ problemId }: { problemId: string }) {
332355 )}
333356 </span>
334357 ))}
335- <label>
336- repeats{" "}
358+ <label
359+ title={
360+ `Each point is timed at least this many times, then repeated ` +
361+ `until it has used ${DEFAULT_TIMING.timeBudgetSeconds} s or ` +
362+ `${DEFAULT_TIMING.maxTimedRuns} runs, and the fastest run is ` +
363+ `reported.`
364+ }
365+ >
366+ min timed runs{" "}
337367 <select
338- value={repeats}
339- onChange={(e) => setRepeats(parseInt(e.target.value, 10))}
368+ value={minRuns}
369+ onChange={(e) => setMinRuns(parseInt(e.target.value, 10))}
340370 >
341- {[1, 3, 5].map((r) => (
371+ {[1, 3, 5, 10].map((r) => (
342372 <option key={r} value={r}>
343373 {r}
344374 </option>
src/app/worker.tsmodified+4−3View file
@@ -7,6 +7,7 @@ import { getSolver } from "../solvers";
77 import { runSweep } from "../harness/sweep";
88 import { runPoint } from "../harness/runner";
99 import { toResultPoint, type ResultPoint } from "../harness/resultSchema";
10+import { DEFAULT_TIMING, type TimingPolicy } from "../harness/timing";
1011 import { matlabBase, solverSource } from "./matlabSources";
1112
1213 export interface SweepRequest {
@@ -14,7 +15,7 @@ export interface SweepRequest {
1415 id: number;
1516 instanceId: string;
1617 solverId: string;
17- repeats: number;
18+ timing: TimingPolicy;
1819 }
1920
2021 export interface SolutionRequest {
@@ -46,7 +47,7 @@ self.onmessage = (e: MessageEvent<WorkerRequest>) => {
4647 instance: getInstance(msg.instanceId),
4748 solver: getSolver(msg.solverId),
4849 sources: { ...matlabBase(), solver: solverSource(msg.solverId) },
49- repeats: msg.repeats,
50+ timing: msg.timing,
5051 onPoint: (p, index, total) => {
5152 const resp: WorkerResponse = {
5253 type: "point",
@@ -68,7 +69,7 @@ self.onmessage = (e: MessageEvent<WorkerRequest>) => {
6869 const p = runPoint({
6970 instance: getInstance(msg.instanceId),
7071 n: msg.n,
71- repeats: 1,
72+ timing: { ...DEFAULT_TIMING, minTimedRuns: 1, timeBudgetSeconds: 0 },
7273 wantGrid: true,
7374 sources: { ...matlabBase(), solver: solverSource(msg.solverId) },
7475 });
src/app/workerClient.tsmodified+3−2View file
@@ -3,6 +3,7 @@
33
44 import type { WorkerRequest, WorkerResponse } from "./worker";
55 import type { ResultPoint } from "../harness/resultSchema";
6+import type { TimingPolicy } from "../harness/timing";
67
78 type Pending = {
89 onPoint?: (point: ResultPoint, index: number, total: number) => void;
@@ -50,13 +51,13 @@ function post(req: WorkerRequest) {
5051 export function sweepInBrowser(
5152 instanceId: string,
5253 solverId: string,
53- repeats: number,
54+ timing: TimingPolicy,
5455 onPoint: (point: ResultPoint, index: number, total: number) => void
5556 ): Promise<ResultPoint[]> {
5657 const id = nextId++;
5758 return new Promise((resolve, reject) => {
5859 pending.set(id, { onPoint, resolve: resolve as never, reject });
59- post({ type: "sweep", id, instanceId, solverId, repeats });
60+ post({ type: "sweep", id, instanceId, solverId, timing });
6061 });
6162 }
6263
src/cli/main.tsmodified+47−13View file
@@ -5,7 +5,8 @@
55 // fastandaccurate list
66 // fastandaccurate run [--instance <id>] [--solver <id>]
77 // [--solver-file f.m --solver-id name]
8-// [--repeats N] [--max-n N] [--label "text"]
8+// [--repeats N] [--time-budget S] [--max-n N]
9+// [--label "text"]
910 // [--out dir]
1011 //
1112 // In development: npx tsx src/cli/main.ts run ...
@@ -17,11 +18,18 @@ import os from "os";
1718 import { INSTANCES, getInstance } from "../problems/laplace2d/spec";
1819 import { setNumblFileIO } from "../harness/numblRun";
1920 import { NodeFileIOAdapter } from "./nodeFileIO";
20-import { ensureChunkie, matlabAvailable, runMatlabSweep } from "./matlabRun";
21+import { matlabAvailable, matlabSetup, runMatlabSweep } from "./matlabRun";
2122
2223 setNumblFileIO((vfs) => new NodeFileIOAdapter(vfs));
23-import { SOLVERS, getSolver, type SolverManifest } from "../solvers";
24+import {
25+ SOLVERS,
26+ getSolver,
27+ solverSourceDir,
28+ sweepNFor,
29+ type SolverManifest,
30+} from "../solvers";
2431 import { runSweep } from "../harness/sweep";
32+import { DEFAULT_TIMING, type TimingPolicy } from "../harness/timing";
2533 import type { MatlabSources } from "../harness/runner";
2634 import {
2735 buildResultFile,
@@ -110,16 +118,29 @@ function listCommand() {
110118 console.log("Problem: laplace-dirichlet-2d (v1)\n");
111119 console.log("Instances:");
112120 for (const inst of INSTANCES) {
113- console.log(` ${inst.id.padEnd(14)} ${inst.label}`);
121+ console.log(` ${inst.id.padEnd(16)} ${inst.label}`);
114122 }
115123 console.log("\nSolvers:");
116124 for (const s of SOLVERS) {
117- console.log(` ${s.id.padEnd(14)} ${s.name} (v${s.version}, ${s.backend})`);
125+ console.log(
126+ ` ${s.id.padEnd(16)} ${s.name} ` +
127+ `(v${s.version}, ${s.backend}, ${s.runtime})`
128+ );
118129 }
119130 }
120131
121132 async function runCommand(flags: Record<string, string>) {
122- const repeats = flags.repeats ? parseInt(flags.repeats, 10) : 5;
133+ const timing: TimingPolicy = {
134+ minTimedRuns: flags.repeats
135+ ? parseInt(flags.repeats, 10)
136+ : DEFAULT_TIMING.minTimedRuns,
137+ timeBudgetSeconds: flags["time-budget"]
138+ ? parseFloat(flags["time-budget"])
139+ : DEFAULT_TIMING.timeBudgetSeconds,
140+ maxTimedRuns: flags["max-repeats"]
141+ ? parseInt(flags["max-repeats"], 10)
142+ : DEFAULT_TIMING.maxTimedRuns,
143+ };
123144 const maxN = flags["max-n"] ? parseInt(flags["max-n"], 10) : undefined;
124145 const outDir = resolve(flags.out ?? "fastandaccurate-results-out");
125146 const instances = flags.instance
@@ -143,7 +164,11 @@ async function runCommand(flags: Record<string, string>) {
143164 version: flags["solver-version"] ?? "0.0.0",
144165 backend: "cpu",
145166 runtime: "numbl",
167+ // A submitted solver sweeps the same resolutions as the reference
168+ // Nystrom solver, per-instance lists included, so its curve lands on
169+ // the same points as the committed ones.
146170 sweepN: getSolver("nystrom-dlp").sweepN,
171+ sweepNByInstance: getSolver("nystrom-dlp").sweepNByInstance,
147172 };
148173 solverList = [
149174 {
@@ -167,7 +192,10 @@ async function runCommand(flags: Record<string, string>) {
167192 }
168193 solverList = wanted.map((manifest) => ({
169194 manifest,
170- sources: { ...base, solver: readSrc(`solvers/${manifest.id}/solver.m`) },
195+ sources: {
196+ ...base,
197+ solver: readSrc(`solvers/${solverSourceDir(manifest)}/solver.m`),
198+ },
171199 source: "builtin",
172200 }));
173201 }
@@ -189,12 +217,14 @@ async function runCommand(flags: Record<string, string>) {
189217 let runEnv = env;
190218 let timer: string | undefined;
191219 if (manifest.runtime === "matlab") {
192- const setup = manifest.id === "chunkie-dlp" ? ensureChunkie() : [];
193- const ns = manifest.sweepN.filter((n) => maxN === undefined || n <= maxN);
220+ const setup = matlabSetup(manifest.id);
221+ const ns = sweepNFor(manifest, inst.id).filter(
222+ (n) => maxN === undefined || n <= maxN
223+ );
194224 const { points, matlabVersion } = runMatlabSweep({
195225 instance: inst,
196226 ns,
197- repeats,
227+ timing,
198228 sources,
199229 setup,
200230 onPoint: printPoint,
@@ -207,7 +237,7 @@ async function runCommand(flags: Record<string, string>) {
207237 instance: inst,
208238 solver: manifest,
209239 sources,
210- repeats,
240+ timing,
211241 maxN,
212242 onPoint: printPoint,
213243 });
@@ -222,7 +252,7 @@ async function runCommand(flags: Record<string, string>) {
222252 source,
223253 },
224254 environment: runEnv,
225- repeats,
255+ timing,
226256 points: resultPoints,
227257 timer,
228258 });
@@ -259,7 +289,11 @@ async function main() {
259289 " --solver-file <f.m> A custom solver file (requires --solver-id)",
260290 " --solver-id <name> Identifier for the custom solver",
261291 " --solver-version <v> Version string for the custom solver",
262- " --repeats <N> Timed repeats per point (default 5)",
292+ " --repeats <N> Minimum timed runs per point (default 5)",
293+ " --time-budget <s> Keep timing a point until it has used this",
294+ " many seconds (default 0.5), which is what",
295+ " makes cheap points reproducible",
296+ " --max-repeats <N> Cap on timed runs per point (default 50)",
263297 " --max-n <N> Restrict the sweep to n <= N",
264298 " --label <text> Free-text machine label recorded in results",
265299 " --out <dir> Output directory (default fastandaccurate-results-out)",
src/cli/matlabRun.tsmodified+32−39View file
@@ -4,18 +4,16 @@
44 // `matlab -batch` once for the whole sweep (one MATLAB startup per
55 // instance), and reads a JSON payload back. Errors are computed on the
66 // node side against the exact solution, as for numbl runs; timing is
7-// MATLAB's own tic/toc with the same warmup-plus-median protocol.
7+// MATLAB's own tic/toc under the same policy as the numbl runner.
88
99 import { execFileSync, spawnSync } from "node:child_process";
10-import { existsSync, mkdirSync, mkdtempSync, readFileSync, rmSync, writeFileSync } from "node:fs";
11-import { homedir, tmpdir } from "node:os";
10+import { existsSync, mkdtempSync, readFileSync, rmSync, writeFileSync } from "node:fs";
11+import { tmpdir } from "node:os";
1212 import { join } from "node:path";
1313 import type { Laplace2dInstance } from "../problems/laplace2d/spec";
1414 import { evalErrors } from "../problems/laplace2d/exact";
1515 import type { ResultPoint } from "../harness/resultSchema";
16-
17-const CHUNKIE_GIT = "https://github.com/fastalgorithms/chunkie";
18-const depsDir = join(homedir(), ".cache", "fastandaccurate", "matlab-deps");
16+import { timedRunLines, type TimingPolicy } from "../harness/timing";
1917
2018 export function matlabAvailable(): boolean {
2119 try {
@@ -26,39 +24,37 @@ export function matlabAvailable(): boolean {
2624 }
2725 }
2826
29-/** Clone chunkie (with the FLAM submodule) on first use and return the
30- * MATLAB setup lines that put it on the path. chunkie's own startup.m is
31- * deliberately not used: it attempts to compile fmm2d when a Fortran
32- * compiler is present, which an unattended run must not do, and the
33- * direct (accel=false) code path needs only the toolbox and FLAM. */
34-export function ensureChunkie(): string[] {
35- const dir = join(depsDir, "chunkie");
36- if (!existsSync(join(dir, "chunkie"))) {
37- mkdirSync(depsDir, { recursive: true });
38- console.log(`fetching chunkie into ${dir}`);
39- execFileSync(
40- "git",
41- ["clone", "--depth", "1", "--recurse-submodules", "--shallow-submodules", CHUNKIE_GIT, dir],
42- { stdio: "inherit" }
43- );
44- }
45- if (!existsSync(join(dir, "chunkie", "FLAM", "startup.m"))) {
46- execFileSync(
47- "git",
48- ["-C", dir, "submodule", "update", "--init", "--depth", "1", "chunkie/FLAM"],
49- { stdio: "inherit" }
50- );
51- }
27+/** The MATLAB setup lines a solver needs before its sweep. Only
28+ * chunkie-dlp has an external dependency; the -mat solvers are plain
29+ * MATLAB and need nothing.
30+ *
31+ * chunkie comes from mip (https://mip.sh), which the harness expects to
32+ * find on the MATLAB path; --install fetches chunkie, and with it FLAM
33+ * and fmm2d, on first use. Taking chunkie from mip rather than from a
34+ * git clone is what makes the accelerated code path available: the mip
35+ * fmm2d package ships a compiled MEX binary for the platform, so
36+ * chunkie's FMM evaluation runs without a Fortran compiler on the
37+ * machine. chunkie's own startup.m is not used; mip puts the three
38+ * packages on the path itself. */
39+export function matlabSetup(solverId: string): string[] {
40+ if (solverId !== "chunkie-dlp") return [];
5241 return [
53- `addpath('${join(dir, "chunkie")}');`,
54- `run('${join(dir, "chunkie", "FLAM", "startup.m")}');`,
42+ // which('/mip') is mip's own idiom for "a non-builtin function named
43+ // mip", so a variable of that name in a user startup cannot mask the
44+ // check.
45+ "if isempty(which('/mip'))",
46+ " error(['chunkie-dlp needs the mip package manager on the MATLAB ' ...",
47+ " 'path. Install it from inside MATLAB with ' ...",
48+ " 'eval(webread(''https://mip.sh/install.txt''))']);",
49+ "end",
50+ "mip load --install chunkie;",
5551 ];
5652 }
5753
5854 export interface MatlabSweepOptions {
5955 instance: Laplace2dInstance;
6056 ns: number[];
61- repeats: number;
57+ timing: TimingPolicy;
6258 sources: { buildProblem: string; bdata: string; solver: string };
6359 /** MATLAB lines run before anything else (addpath etc.). */
6460 setup: string[];
@@ -75,7 +71,7 @@ function asArray(x: number | number[]): number[] {
7571 }
7672
7773 export function runMatlabSweep(opts: MatlabSweepOptions): MatlabSweepResult {
78- const { instance, ns, repeats } = opts;
74+ const { instance, ns, timing } = opts;
7975 const dir = mkdtempSync(join(tmpdir(), "fastandaccurate-matlab-"));
8076 try {
8177 writeFileSync(join(dir, "build_problem.m"), opts.sources.buildProblem);
@@ -85,8 +81,8 @@ export function runMatlabSweep(opts: MatlabSweepOptions): MatlabSweepResult {
8581 "% generated by the fastandaccurate MATLAB harness",
8682 ...opts.setup,
8783 `ns = [${ns.join(" ")}];`,
88- `nrep = ${repeats};`,
89- `prob = build_problem(${instance.a}, ${instance.k}, ${instance.d}, 0);`,
84+ `prob = build_problem(${instance.a}, ${instance.k}, ${instance.p ?? 0}, ` +
85+ `${instance.d}, 0, ${instance.nearBoundary ? 1 : 0});`,
9086 "results = cell(numel(ns), 1);",
9187 "% Session-level warmup: the whole sweep shares one MATLAB process, so",
9288 "% without this the first resolution absorbs all of the one-time cost",
@@ -98,10 +94,7 @@ export function runMatlabSweep(opts: MatlabSweepOptions): MatlabSweepResult {
9894 " n = ns(i);",
9995 " tic; out = solver(prob, n); cold = toc;",
10096 " out = solver(prob, n);",
101- " times = zeros(nrep, 1);",
102- " for r = 1:nrep",
103- " tic; out = solver(prob, n); times(r) = toc;",
104- " end",
97+ ...timedRunLines("out = solver(prob, n)", "times", timing).map((l) => ` ${l}`),
10598 " results{i} = struct('n', n, 'cold', cold, 'times', times, 'ueval', out.uEval);",
10699 " fprintf('point n=%d done (%.3fs)\\n', n, min(times));",
107100 "end",
src/harness/resultSchema.tsmodified+15−3View file
@@ -10,6 +10,7 @@ import {
1010 type Laplace2dInstance,
1111 } from "../problems/laplace2d/spec";
1212 import type { RunPoint } from "./runner";
13+import type { TimingPolicy } from "./timing";
1314
1415 export const RESULT_FORMAT = "fastandaccurate-result";
1516 export const RESULT_FORMAT_VERSION = 1;
@@ -57,10 +58,19 @@ export interface ResultFile {
5758 environment: ResultEnvironment;
5859 protocol: {
5960 warmupRuns: number;
60- timedRuns: number;
61+ /** The adaptive timed-run policy: at least minTimedRuns runs, then as
62+ * many more as fit in timeBudgetSeconds, never past maxTimedRuns. How
63+ * many a given point actually ran is the length of its
64+ * solveSecondsAll. */
65+ minTimedRuns: number;
66+ timeBudgetSeconds: number;
67+ maxTimedRuns: number;
6168 /** Which statistic of the timed runs is reported as solveSeconds. */
6269 statistic: "min";
6370 timer: string;
71+ /** Results written before the adaptive policy carry a fixed count
72+ * here instead of the three fields above. */
73+ timedRuns?: number;
6474 };
6575 createdUtc: string;
6676 points: ResultPoint[];
@@ -92,7 +102,7 @@ export async function buildResultFile(opts: {
92102 instance: Laplace2dInstance;
93103 solver: { id: string; version: string; backend: "cpu" | "gpu"; source: string };
94104 environment: ResultEnvironment;
95- repeats: number;
105+ timing: TimingPolicy;
96106 points: ResultPoint[];
97107 /** What measured the times (default numbl tic/toc). */
98108 timer?: string;
@@ -109,7 +119,9 @@ export async function buildResultFile(opts: {
109119 environment: opts.environment,
110120 protocol: {
111121 warmupRuns: 2,
112- timedRuns: opts.repeats,
122+ minTimedRuns: opts.timing.minTimedRuns,
123+ timeBudgetSeconds: opts.timing.timeBudgetSeconds,
124+ maxTimedRuns: opts.timing.maxTimedRuns,
113125 statistic: "min",
114126 timer: opts.timer ?? "numbl tic/toc",
115127 },
src/harness/runner.tsmodified+18−15View file
@@ -1,14 +1,17 @@
11 // Runs one solver at one resolution on one instance and measures it.
22 // The timing protocol (see docs/problems/laplace-dirichlet-2d.md): two
3-// untimed warmup runs absorb JIT compilation, then `repeats` timed runs
4-// of which the fastest is reported. The minimum is the least contaminated
5-// estimator of the solver's own cost, since everything that interferes
6-// (scheduling, residual compilation, other load) only ever adds time; all
7-// the individual timings are recorded too. The first warmup is timed and
8-// reported separately as the cold time. All timing is MATLAB tic/toc
9-// inside the numbl session, so browser and node measure the same thing.
3+// untimed warmup runs absorb JIT compilation, then timed runs under the
4+// policy in ./timing (a floor on the count, then as many more as fit in a
5+// time budget), of which the fastest is reported. The minimum is the least
6+// contaminated estimator of the solver's own cost, since everything that
7+// interferes (scheduling, residual compilation, other load) only ever adds
8+// time; all the individual timings are recorded too. The first warmup is
9+// timed and reported separately as the cold time. All timing is MATLAB
10+// tic/toc inside the numbl session, so browser and node measure the same
11+// thing.
1012
1113 import { runNumblScript } from "./numblRun";
14+import { DEFAULT_TIMING, timedRunLines, type TimingPolicy } from "./timing";
1215 import type { Laplace2dInstance } from "../problems/laplace2d/spec";
1316 import { evalErrors } from "../problems/laplace2d/exact";
1417
@@ -24,8 +27,8 @@ export interface MatlabSources {
2427 export interface RunPointRequest {
2528 instance: Laplace2dInstance;
2629 n: number;
27- /** Timed repeats after the warmups (default 5). */
28- repeats?: number;
30+ /** How the timed runs are counted (default DEFAULT_TIMING). */
31+ timing?: TimingPolicy;
2932 /** Also evaluate the solution on the visualization grid. */
3033 wantGrid?: boolean;
3134 sources: MatlabSources;
@@ -35,6 +38,7 @@ export interface RunPoint {
3538 n: number;
3639 /** The fastest of the timed runs. */
3740 solveSeconds: number;
41+ /** Every timed run, in order; its length is how many the policy ran. */
3842 solveSecondsAll: number[];
3943 coldSeconds: number;
4044 relMax: number;
@@ -50,17 +54,16 @@ function numLiteral(x: number): string {
5054
5155 export function runPoint(req: RunPointRequest): RunPoint {
5256 const { instance, n } = req;
53- const repeats = req.repeats ?? 5;
57+ const timing = req.timing ?? DEFAULT_TIMING;
5458 const wantGrid = req.wantGrid ?? false;
5559 const main = [
5660 "% generated by the fastandaccurate harness",
57- `prob = build_problem(${numLiteral(instance.a)}, ${numLiteral(instance.k)}, ${numLiteral(instance.d)}, ${wantGrid ? 1 : 0});`,
61+ `prob = build_problem(${numLiteral(instance.a)}, ${numLiteral(instance.k)}, ` +
62+ `${numLiteral(instance.p ?? 0)}, ${numLiteral(instance.d)}, ${wantGrid ? 1 : 0}, ` +
63+ `${instance.nearBoundary ? 1 : 0});`,
5864 `tic; out = solver(prob, ${n}); res_cold = toc;`,
5965 `out = solver(prob, ${n});`,
60- `res_times = zeros(${repeats}, 1);`,
61- `for irep = 1:${repeats}`,
62- ` tic; out = solver(prob, ${n}); res_times(irep) = toc;`,
63- "end",
66+ ...timedRunLines(`out = solver(prob, ${n})`, "res_times", timing),
6467 "res_ueval = out.uEval;",
6568 "res_ugrid = out.uGrid;",
6669 "",
src/harness/sweep.tsmodified+5−4View file
@@ -2,21 +2,22 @@
22 // one instance, reporting each point as it lands.
33
44 import { runPoint, type MatlabSources, type RunPoint } from "./runner";
5+import { DEFAULT_TIMING, type TimingPolicy } from "./timing";
56 import type { Laplace2dInstance } from "../problems/laplace2d/spec";
6-import type { SolverManifest } from "../solvers";
7+import { sweepNFor, type SolverManifest } from "../solvers";
78
89 export interface SweepOptions {
910 instance: Laplace2dInstance;
1011 solver: SolverManifest;
1112 sources: MatlabSources;
12- repeats?: number;
13+ timing?: TimingPolicy;
1314 /** Restrict the sweep to n values <= this (for quick runs). */
1415 maxN?: number;
1516 onPoint?: (point: RunPoint, index: number, total: number) => void;
1617 }
1718
1819 export function runSweep(opts: SweepOptions): RunPoint[] {
19- const ns = opts.solver.sweepN.filter(
20+ const ns = sweepNFor(opts.solver, opts.instance.id).filter(
2021 (n) => opts.maxN === undefined || n <= opts.maxN
2122 );
2223 const points: RunPoint[] = [];
@@ -24,7 +25,7 @@ export function runSweep(opts: SweepOptions): RunPoint[] {
2425 const p = runPoint({
2526 instance: opts.instance,
2627 n,
27- repeats: opts.repeats ?? 5,
28+ timing: opts.timing ?? DEFAULT_TIMING,
2829 sources: opts.sources,
2930 });
3031 points.push(p);
src/harness/timing.tsadded+55−0View file
@@ -0,0 +1,55 @@
1+// The timing policy, shared by the numbl runner and the MATLAB runner so
2+// that both measure the same way (see docs/problems/laplace-dirichlet-2d.md).
3+//
4+// A fixed repeat count times a cheap point badly: a solve of 0.1 ms is
5+// dominated by whatever the operating system was doing during those few
6+// hundred microseconds, and the min of five such samples still scatters by
7+// a factor of several from one sweep to the next, which is what made the
8+// low-resolution end of the work-precision curves noisy. So the count is
9+// adaptive: keep running timed repetitions until they have accumulated a
10+// real amount of work, subject to a floor on the count and a cap so that
11+// an expensive solve is not repeated forever. The reported time is still
12+// the minimum over the timed runs, and every timing is recorded.
13+
14+export interface TimingPolicy {
15+ /** Timed runs never stop below this count. */
16+ minTimedRuns: number;
17+ /** Timed runs continue past the minimum until they have accumulated
18+ * this much time in total. */
19+ timeBudgetSeconds: number;
20+ /** Hard cap on the count. It bounds both the time a cheap point can
21+ * spend and the number of timings a result file has to carry. */
22+ maxTimedRuns: number;
23+}
24+
25+export const DEFAULT_TIMING: TimingPolicy = {
26+ minTimedRuns: 5,
27+ timeBudgetSeconds: 0.5,
28+ maxTimedRuns: 50,
29+};
30+
31+/**
32+ * MATLAB lines that run `call` under the policy and leave one timing per
33+ * run in `timesVar` (a column vector, trimmed to the number actually
34+ * run). The loop lives on the MATLAB side so that all timing is tic/toc
35+ * inside the solver's own runtime.
36+ */
37+export function timedRunLines(
38+ call: string,
39+ timesVar: string,
40+ policy: TimingPolicy
41+): string[] {
42+ const { minTimedRuns: min, maxTimedRuns: max, timeBudgetSeconds: budget } = policy;
43+ return [
44+ `${timesVar} = zeros(${max}, 1);`,
45+ "res_nrun = 0;",
46+ "res_tot = 0;",
47+ `while res_nrun < ${max} && (res_nrun < ${min} || res_tot < ${budget})`,
48+ ` tic; ${call}; res_one = toc;`,
49+ " res_nrun = res_nrun + 1;",
50+ ` ${timesVar}(res_nrun) = res_one;`,
51+ " res_tot = res_tot + res_one;",
52+ "end",
53+ `${timesVar} = ${timesVar}(1:res_nrun);`,
54+ ];
55+}
src/problems/laplace2d/exact.tsmodified+93−10View file
@@ -3,7 +3,7 @@
33 // implementations of these simple formulas check each other: a solver can
44 // only reach high accuracy if both agree.
55
6-import type { Laplace2dInstance } from "./spec";
6+import { hasCorners, hasNearBoundary, type Laplace2dInstance } from "./spec";
77
88 export interface Source {
99 x: number;
@@ -11,26 +11,93 @@ export interface Source {
1111 c: number;
1212 }
1313
14-/** Boundary radius r(t) = 1 + a cos(k t). */
14+/** Boundary radius: r(t) = 1 + a cos(k t) for the star family, and
15+ * r(t) = (cos^p t + sin^p t)^(-1/p) for the rounded square, which traces
16+ * the superellipse |x|^p + |y|^p = 1. Written through a logarithm because
17+ * cos^p t underflows for p in the hundreds. */
1518 export function boundaryR(inst: Laplace2dInstance, t: number): number {
19+ if (inst.shape === "rounded-square") {
20+ const p = inst.p as number;
21+ return Math.exp(-Math.log(Math.cos(t) ** p + Math.sin(t) ** p) / p);
22+ }
1623 return 1 + inst.a * Math.cos(inst.k * t);
1724 }
1825
26+/** dr/dt of the boundary radius. Only the ratio f'/f enters, which keeps
27+ * the rounded-square case away from the underflow in f itself. */
28+export function boundaryRD(inst: Laplace2dInstance, t: number): number {
29+ if (inst.shape === "rounded-square") {
30+ const p = inst.p as number;
31+ const c = Math.cos(t);
32+ const sn = Math.sin(t);
33+ const f = c ** p + sn ** p;
34+ const fp = p * (sn ** (p - 1) * c - c ** (p - 1) * sn);
35+ return (-boundaryR(inst, t) / p) * (fp / f);
36+ }
37+ return -inst.a * inst.k * Math.sin(inst.k * t);
38+}
39+
40+/** The largest radius the boundary reaches, which sets the view extent
41+ * and the visualization grid. */
42+export function maxRadius(inst: Laplace2dInstance): number {
43+ if (inst.shape === "rounded-square") return boundaryR(inst, Math.PI / 4);
44+ return 1 + Math.abs(inst.a);
45+}
46+
47+/** Angles of the corners, for instances that have them: the four
48+ * diagonals of the rounded square. */
49+export function cornerAngles(inst: Laplace2dInstance): number[] {
50+ if (!hasCorners(inst)) return [];
51+ return [0, 1, 2, 3].map((j) => Math.PI / 4 + (j * Math.PI) / 2);
52+}
53+
54+/** Distances inside the boundary of the near-field evaluation points. The
55+ * same four distances serve the near-corner set of an instance with
56+ * corners and the near-boundary set of an instance that carries one: the
57+ * smallest is a third of the corner radius at p = 100, and about half a
58+ * node spacing at the resolutions the sweeps reach. */
59+export const NEAR_DELTAS = [0.005, 0.01, 0.02, 0.05];
60+
61+/** Boundary parameters of the near-boundary evaluation points: eight
62+ * around the curve, offset so that they fall at no special phase of a lobe
63+ * pattern. Empty on an instance without the set. */
64+export function nearBoundaryParams(inst: Laplace2dInstance): number[] {
65+ if (!hasNearBoundary(inst)) return [];
66+ return Array.from({ length: 8 }, (_, j) => (2 * Math.PI * j) / 8 + 0.07);
67+}
68+
1969 /** Boundary point at parameter t. */
2070 export function boundaryPoint(inst: Laplace2dInstance, t: number) {
2171 const r = boundaryR(inst, t);
2272 return { x: r * Math.cos(t), y: r * Math.sin(t) };
2373 }
2474
75+/** The point a distance delta inside the boundary along the inward unit
76+ * normal at parameter t. Since delta stays below the smallest radius of
77+ * curvature on the instances that use this, the distance from the point to
78+ * the curve is exactly delta. */
79+export function inwardPoint(inst: Laplace2dInstance, t: number, delta: number) {
80+ const r = boundaryR(inst, t);
81+ const dr = boundaryRD(inst, t);
82+ const dx = dr * Math.cos(t) - r * Math.sin(t);
83+ const dy = dr * Math.sin(t) + r * Math.cos(t);
84+ const sp = Math.hypot(dx, dy);
85+ // The outward unit normal of the counterclockwise curve is (y', -x')/|x'|.
86+ return {
87+ x: r * Math.cos(t) - (delta * dy) / sp,
88+ y: r * Math.sin(t) + (delta * dx) / sp,
89+ };
90+}
91+
2592 /** The three exact-solution sources: boundary points at phi_j pushed a
2693 * distance d along the outward normal, strengths [1.0, -0.6, 0.8]. */
2794 export function sources(inst: Laplace2dInstance): Source[] {
28- const { a, k, d } = inst;
95+ const { d } = inst;
2996 const strengths = [1.0, -0.6, 0.8];
3097 return strengths.map((c, j) => {
3198 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);
99+ const r = boundaryR(inst, phi);
100+ const dr = boundaryRD(inst, phi);
34101 const bx = r * Math.cos(phi);
35102 const by = r * Math.sin(phi);
36103 const dx = dr * Math.cos(phi) - r * Math.sin(phi);
@@ -49,9 +116,14 @@ export function exactU(inst: Laplace2dInstance, x: number, y: number): number {
49116 return u;
50117 }
51118
52-/** The 289 evaluation points: 32 rays, radial fractions 0.1..0.9, plus
53- * the origin. Order matches build_problem.m: radius outer, angle inner,
54- * origin last. */
119+/** The evaluation points: 32 rays, radial fractions 0.1..0.9, plus the
120+ * origin (289 points); on an instance with corners a further four points
121+ * per corner just inside the boundary along its diagonal (305 in all); and
122+ * on an instance carrying the near-boundary set four points per parameter
123+ * along the inward normal (321 in all). Order matches build_problem.m:
124+ * radius outer, angle inner, origin last, then the near-corner points with
125+ * the corner index outer and the distance inner, then the near-boundary
126+ * points with the parameter outer and the distance inner. */
55127 export function evalPoints(inst: Laplace2dInstance): { x: number; y: number }[] {
56128 const rho = [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9];
57129 const pts: { x: number; y: number }[] = [];
@@ -63,6 +135,17 @@ export function evalPoints(inst: Laplace2dInstance): { x: number; y: number }[]
63135 }
64136 }
65137 pts.push({ x: 0, y: 0 });
138+ for (const th of cornerAngles(inst)) {
139+ const rc = boundaryR(inst, th);
140+ for (const delta of NEAR_DELTAS) {
141+ pts.push({ x: (rc - delta) * Math.cos(th), y: (rc - delta) * Math.sin(th) });
142+ }
143+ }
144+ for (const t of nearBoundaryParams(inst)) {
145+ for (const delta of NEAR_DELTAS) {
146+ pts.push(inwardPoint(inst, t, delta));
147+ }
148+ }
66149 return pts;
67150 }
68151
@@ -75,13 +158,13 @@ export function exactAtEvalPoints(inst: Laplace2dInstance): Float64Array {
75158 }
76159
77160 /** 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],
161+ * R = 1.05 max_t r(t). Flat index p = ix * ngrid + iy with x = xs[ix],
79162 * y = xs[iy] (y varies fastest), matching build_problem.m's meshgrid
80163 * column order. */
81164 export const VIZ_NGRID = 200;
82165
83166 export function vizGrid(inst: Laplace2dInstance) {
84- const R = 1.05 * (1 + Math.abs(inst.a));
167+ const R = 1.05 * maxRadius(inst);
85168 const xs = new Float64Array(VIZ_NGRID);
86169 for (let i = 0; i < VIZ_NGRID; i++) {
87170 xs[i] = -R + (2 * R * i) / (VIZ_NGRID - 1);
src/problems/laplace2d/matlab/build_problem.mmodified+149−20View file
@@ -1,10 +1,24 @@
1-function prob = build_problem(a, k, d, wantGrid)
1+function prob = build_problem(a, k, p, d, wantGrid, nearBoundary)
22 % BUILD_PROBLEM Assemble the problem struct for laplace-dirichlet-2d.
33 %
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).
4+% The domain is bounded by one of two analytic curves, selected by p:
5+%
6+% p == 0 the star family, r(t) = 1 + a cos(k t),
7+% p > 0 the rounded square, r(t) = (cos^p t + sin^p t)^(-1/p), which
8+% traces the superellipse |x|^p + |y|^p = 1: four nearly
9+% straight sides meeting at four corners of radius about 1.4/p.
10+% a and k are unused. p must be even, which keeps the curve
11+% analytic, so curveDD exists everywhere as the interface
12+% requires.
13+%
14+% In both cases x(t) = r(t) [cos t; sin t], t in [0, 2 pi).
15+%
16+% nearBoundary (default 0) adds the near-boundary target set to the
17+% evaluation points: four points per parameter just inside the boundary
18+% along the inward normal, at eight parameters around the curve.
19+%
620 % 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
21+% sources outside the domain: s_j is the boundary point at parameter
822 % phi_j = 2 pi (j-1)/3 + 0.4 pushed a distance d along the outward unit
923 % normal, with strengths c = [1.0; -0.6; 0.8].
1024 %
@@ -18,44 +32,111 @@ function prob = build_problem(a, k, d, wantGrid)
1832 % curveD @(t) -> [x' y'] first derivative
1933 % curveDD @(t) -> [x'' y''] second derivative
2034 % g @(t) -> g Dirichlet data at boundary parameter t
21-% evalXY 289 x 2 points where uEval is required
35+% evalXY m x 2 points where uEval is required (289; 305
36+% on an instance with corners, 321 on one
37+% with the near-boundary set)
2238 % vizXY m x 2 grid points where uGrid is requested
2339 % (m = 0 when no visualization is wanted)
2440
41+if nargin < 6
42+ nearBoundary = 0;
43+end
44+
45+corners = p > 0;
46+
47+if corners
48+ prob = struct();
49+ prob.curve = @(t) sq_xy(t, p);
50+ prob.curveD = @(t) sq_xyd(t, p);
51+ prob.curveDD = @(t) sq_xydd(t, p);
52+ radius = @(t) sq_r(t, p);
53+ radiusD = @(t) sq_rd(t, p);
54+ rmax = sq_r(pi/4, p);
55+else
56+ prob = struct();
57+ prob.curve = @(t) [(1 + a*cos(k*t)).*cos(t), (1 + a*cos(k*t)).*sin(t)];
58+ prob.curveD = @(t) [-a*k*sin(k*t).*cos(t) - (1 + a*cos(k*t)).*sin(t), ...
59+ -a*k*sin(k*t).*sin(t) + (1 + a*cos(k*t)).*cos(t)];
60+ prob.curveDD = @(t) [(-a*k*k*cos(k*t) - 1 - a*cos(k*t)).*cos(t) + 2*a*k*sin(k*t).*sin(t), ...
61+ (-a*k*k*cos(k*t) - 1 - a*cos(k*t)).*sin(t) - 2*a*k*sin(k*t).*cos(t)];
62+ radius = @(t) 1 + a*cos(k*t);
63+ radiusD = @(t) -a*k*sin(k*t);
64+ rmax = 1 + abs(a);
65+end
66+
67+% The three sources: boundary points at phi_j pushed d along the outward
68+% normal. The normal comes from the radius and its derivative, so this is
69+% the same formula for both families.
2570 phi = 2*pi*[0; 1; 2]/3 + 0.4;
2671 c = [1.0; -0.6; 0.8];
27-rphi = 1 + a*cos(k*phi);
72+rphi = radius(phi);
73+drphi = radiusD(phi);
2874 bx = rphi.*cos(phi);
2975 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);
76+dxb = drphi.*cos(phi) - rphi.*sin(phi);
77+dyb = drphi.*sin(phi) + rphi.*cos(phi);
3278 sp = sqrt(dxb.^2 + dyb.^2);
3379 sx = bx + d*(dyb./sp);
3480 sy = by - d*(dxb./sp);
3581
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);
82+prob.g = @(t) laplace2d_bdata(t, radius(t), sx, sy, c);
4383
4484 % Evaluation points: 32 rays, radial fractions 0.1..0.9, plus the origin
45-% (289 points). The rule must match evalPoints() in
46-% src/problems/laplace2d/exact.ts.
85+% (289 points). On an instance with corners, four further points per
86+% corner follow, just inside the boundary along the corner's diagonal
87+% (305 points). On an instance carrying the near-boundary set, four
88+% further points per parameter follow after those, at the same four
89+% distances inside the boundary but along the inward unit normal, at
90+% eight parameters around the curve (321 points). The rule must match
91+% evalPoints() in src/problems/laplace2d/exact.ts.
4792 rho = (1:9)'/10;
4893 th = 2*pi*(0:31)'/32 + 0.13;
49-pts = zeros(numel(rho)*numel(th) + 1, 2);
94+deltas = [0.005; 0.01; 0.02; 0.05];
95+if corners
96+ cth = pi/4 + pi*(0:3)'/2;
97+else
98+ cth = zeros(0, 1);
99+end
100+if nearBoundary
101+ ntt = 2*pi*(0:7)'/8 + 0.07;
102+else
103+ ntt = zeros(0, 1);
104+end
105+npts = numel(rho)*numel(th) + 1 + (numel(cth) + numel(ntt))*numel(deltas);
106+pts = zeros(npts, 2);
50107 idx = 1;
51108 for i = 1:numel(rho)
52109 for j = 1:numel(th)
53- rr = rho(i)*(1 + a*cos(k*th(j)));
110+ rr = rho(i)*radius(th(j));
54111 pts(idx, 1) = rr*cos(th(j));
55112 pts(idx, 2) = rr*sin(th(j));
56113 idx = idx + 1;
57114 end
58115 end
116+pts(idx, 1) = 0;
117+pts(idx, 2) = 0;
118+idx = idx + 1;
119+for i = 1:numel(cth)
120+ rc = radius(cth(i));
121+ for j = 1:numel(deltas)
122+ pts(idx, 1) = (rc - deltas(j))*cos(cth(i));
123+ pts(idx, 2) = (rc - deltas(j))*sin(cth(i));
124+ idx = idx + 1;
125+ end
126+end
127+% The inward normal of the counterclockwise curve is -(y', -x')/|x'|.
128+rn = radius(ntt);
129+drn = radiusD(ntt);
130+dxn = drn.*cos(ntt) - rn.*sin(ntt);
131+dyn = drn.*sin(ntt) + rn.*cos(ntt);
132+spn = sqrt(dxn.^2 + dyn.^2);
133+for i = 1:numel(ntt)
134+ for j = 1:numel(deltas)
135+ pts(idx, 1) = rn(i)*cos(ntt(i)) - deltas(j)*dyn(i)/spn(i);
136+ pts(idx, 2) = rn(i)*sin(ntt(i)) + deltas(j)*dxn(i)/spn(i);
137+ idx = idx + 1;
138+ end
139+end
59140 prob.evalXY = pts;
60141
61142 % Visualization grid: ngrid x ngrid points over the bounding square,
@@ -63,7 +144,7 @@ prob.evalXY = pts;
63144 % the domain are included; the viewer masks them.
64145 if wantGrid
65146 ngrid = 200;
66- R = 1.05*(1 + abs(a));
147+ R = 1.05*rmax;
67148 xs = linspace(-R, R, ngrid);
68149 [X, Y] = meshgrid(xs, xs);
69150 prob.vizXY = [X(:), Y(:)];
@@ -72,3 +153,51 @@ else
72153 end
73154
74155 end
156+
157+% --- the rounded-square radius and its derivatives ---------------------
158+%
159+% With f(t) = cos^p t + sin^p t, the radius is r = f^(-1/p). For p in the
160+% hundreds f underflows near the diagonals, so r goes through a logarithm
161+% and the derivatives are written with the ratios f'/f and f''/f, which
162+% stay of moderate size.
163+
164+function r = sq_r(t, p)
165+c = cos(t); s = sin(t);
166+r = exp(-log(c.^p + s.^p)/p);
167+end
168+
169+function r1 = sq_rd(t, p)
170+c = cos(t); s = sin(t);
171+f = c.^p + s.^p;
172+fp = p*(s.^(p-1).*c - c.^(p-1).*s);
173+r1 = -(sq_r(t, p)/p).*(fp./f);
174+end
175+
176+function r2 = sq_rdd(t, p)
177+c = cos(t); s = sin(t);
178+f = c.^p + s.^p;
179+fp = p*(s.^(p-1).*c - c.^(p-1).*s);
180+fpp = p*((p-1)*(s.^(p-2).*c.^2 + c.^(p-2).*s.^2) - f);
181+u = fp./f;
182+v = fpp./f;
183+r2 = sq_r(t, p).*((1/p)*(1/p + 1)*u.^2 - v/p);
184+end
185+
186+function xy = sq_xy(t, p)
187+r = sq_r(t, p);
188+xy = [r.*cos(t), r.*sin(t)];
189+end
190+
191+function xy = sq_xyd(t, p)
192+r = sq_r(t, p);
193+r1 = sq_rd(t, p);
194+xy = [r1.*cos(t) - r.*sin(t), r1.*sin(t) + r.*cos(t)];
195+end
196+
197+function xy = sq_xydd(t, p)
198+r = sq_r(t, p);
199+r1 = sq_rd(t, p);
200+r2 = sq_rdd(t, p);
201+xy = [r2.*cos(t) - 2*r1.*sin(t) - r.*cos(t), ...
202+ r2.*sin(t) + 2*r1.*cos(t) - r.*sin(t)];
203+end
src/problems/laplace2d/matlab/laplace2d_bdata.mmodified+5−5View file
@@ -1,9 +1,9 @@
1-function g = laplace2d_bdata(t, a, k, sx, sy, c)
1+function g = laplace2d_bdata(t, r, sx, sy, c)
22 % 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);
3+% Evaluates u(x(t)) = sum_j c_j log|x(t) - s_j| at boundary parameters t,
4+% given the boundary radius r = r(t) at those parameters (which family the
5+% radius came from does not matter here). Called through the prob.g handle
6+% built in build_problem; solvers see only that handle.
77 x = r.*cos(t);
88 y = r.*sin(t);
99 g = zeros(size(t));
src/problems/laplace2d/spec.tsmodified+95−2View file
@@ -11,14 +11,39 @@ export interface Laplace2dInstance {
1111 /** Short stable identifier used in results and URLs. */
1212 id: string;
1313 label: string;
14- /** Boundary r(t) = 1 + a cos(k t). */
14+ /** Which boundary family the instance draws from. "star", the default,
15+ * is r(t) = 1 + a cos(k t). "rounded-square" is the superellipse
16+ * |x|^p + |y|^p = 1, four nearly straight sides meeting at four rounded
17+ * corners whose radius falls off like 1.4/p; it is analytic for even p,
18+ * so the interface a solver sees is unchanged. */
19+ shape?: "star" | "rounded-square";
20+ /** Boundary r(t) = 1 + a cos(k t) for the star family; both are 0 and
21+ * unused for rounded-square. */
1522 a: number;
1623 k: number;
24+ /** Corner sharpness of the rounded-square family (an even integer). */
25+ p?: number;
1726 /** Distance of the exact solution's sources beyond the boundary. */
1827 d: number;
28+ /** Whether the instance carries the near-boundary target set in addition
29+ * to the standard evaluation points: points a few thousandths inside the
30+ * boundary along the inward normal, where a quadrature rule with no
31+ * near-field correction loses its accuracy. */
32+ nearBoundary?: boolean;
1933 description: string;
2034 }
2135
36+/** Whether the instance has corners, and so carries near-corner
37+ * evaluation points in addition to the standard set. */
38+export function hasCorners(inst: Laplace2dInstance): boolean {
39+ return inst.shape === "rounded-square";
40+}
41+
42+/** Whether the instance carries the near-boundary target set. */
43+export function hasNearBoundary(inst: Laplace2dInstance): boolean {
44+ return inst.nearBoundary === true;
45+}
46+
2247 export const INSTANCES: Laplace2dInstance[] = [
2348 {
2449 id: "disk-easy",
@@ -52,6 +77,64 @@ export const INSTANCES: Laplace2dInstance[] = [
5277 "data barely continues past the boundary, which defeats methods " +
5378 "whose representation assumes it does.",
5479 },
80+ {
81+ id: "flower-15",
82+ label: "15-lobe flower, complex geometry",
83+ a: 0.2,
84+ k: 15,
85+ d: 0.25,
86+ description:
87+ "Fifteen lobes, with boundary curvature up to 69 (radius 0.014) " +
88+ "against 13.9 on star-hard. Here the sources sit a full 0.25 from " +
89+ "the boundary, so the data continues comfortably and the difficulty " +
90+ "is the geometry alone: a method must resolve a boundary whose " +
91+ "features are twenty times smaller than the domain, and the " +
92+ "evaluation points, at a fixed fraction of the local radius, fall " +
93+ "within 0.03 of it because the boundary runs nearly radially along " +
94+ "the flanks of the lobes.",
95+ },
96+ {
97+ id: "square-corners",
98+ label: "Rounded square, targets in the corners",
99+ shape: "rounded-square",
100+ p: 100,
101+ a: 0,
102+ k: 0,
103+ d: 0.25,
104+ description:
105+ "A square with corners rounded to a radius of 0.014, one percent of " +
106+ "its side, and sides that are straight to within a thousandth. " +
107+ "Unlike the star instances, the difficulty is local: four small " +
108+ "patches of the boundary need everything the discretization has, and " +
109+ "the rest needs almost nothing. This instance also carries sixteen " +
110+ "extra evaluation points, four at each corner, at 0.005, 0.01, 0.02 " +
111+ "and 0.05 inside along the diagonal. The closest sits well within " +
112+ "the corner's own radius, which is where a quadrature rule with no " +
113+ "near-field correction gives up: the reported error on this instance " +
114+ "is normally attained at those points rather than in the bulk.",
115+ },
116+ {
117+ id: "star-nearfield",
118+ label: "5-lobe star, targets against the boundary",
119+ a: 0.3,
120+ k: 5,
121+ d: 0.08,
122+ nearBoundary: true,
123+ description:
124+ "The domain and the data of star-hard, asked a different question. " +
125+ "Besides the standard evaluation points it carries 32 targets 0.005 " +
126+ "to 0.05 inside the boundary, along the inward normal at eight " +
127+ "places around the curve. The closest of them is about half a node " +
128+ "spacing from the boundary at the resolutions these sweeps reach, " +
129+ "and that is where a quadrature rule with no near-field correction " +
130+ "stops converging: the error of the plain periodic trapezoid rule " +
131+ "at a target a distance delta inside behaves like " +
132+ "exp(-2 pi delta / h) in the node spacing h, so it can hold in the " +
133+ "bulk long after it has failed here. Because the geometry and the " +
134+ "sources are exactly those of star-hard, the difference between the " +
135+ "two instances measures one thing only: whether a method can " +
136+ "evaluate its own representation close to the boundary.",
137+ },
55138 ];
56139
57140 /** The instance a visitor sees first: the one that separates the methods
@@ -70,7 +153,7 @@ export function getInstance(id: string): Laplace2dInstance {
70153 * future artifact cache hashes it).
71154 */
72155 export function canonicalSpec(inst: Laplace2dInstance) {
73- return {
156+ const base = {
74157 a: inst.a,
75158 d: inst.d,
76159 instance: inst.id,
@@ -78,6 +161,16 @@ export function canonicalSpec(inst: Laplace2dInstance) {
78161 problem: PROBLEM_ID,
79162 problemVersion: PROBLEM_VERSION,
80163 };
164+ // The shape keys appear only for instances that are not from the star
165+ // family, and nearBoundary only where it is set, so that the spec
166+ // strings and hashes of the instances that existed before each of those
167+ // fields are exactly what they were.
168+ const shaped =
169+ !inst.shape || inst.shape === "star"
170+ ? base
171+ : { ...base, p: inst.p, shape: inst.shape };
172+ if (!inst.nearBoundary) return shaped;
173+ return { ...shaped, nearBoundary: true };
81174 }
82175
83176 export function canonicalSpecJson(inst: Laplace2dInstance): string {
src/solvers/chunkie-dlp/solver.mmodified+8−15View file
@@ -11,8 +11,8 @@ function out = solver(prob, n)
1111 % with corrected quadrature for targets near the boundary.
1212 %
1313 % This solver runs in real MATLAB only: the command-line harness invokes
14-% it through `matlab -batch` with chunkie on the path (fetched on first
15-% use). It is not runnable in the browser.
14+% it through `matlab -batch`, having put chunkie on the path with
15+% `mip load --install chunkie`. It is not runnable in the browser.
1616 %
1717 % n : number of chunks (16 points each).
1818
@@ -41,19 +41,12 @@ end
4141 end
4242
4343 function u = eval_targets(chnkr, fkern, sigma, XY)
44-% Direct (unaccelerated) evaluation, in blocks to bound memory. accel is
45-% disabled so that the FLAM/fmm2d submodules are not required; at these
46-% 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
44+% Evaluation with chunkie's own defaults, which means its FMM
45+% acceleration of the smooth part (opts.accel, true unless overridden)
46+% wherever the point counts make that worth doing. The mip fmm2d package
47+% supplies the compiled binary this needs.
48+u = chunkerkerneval(chnkr, fkern, sigma, XY.');
49+u = u(:);
5750 end
5851
5952 function [r, d, d2] = fcurve(t, prob)
src/solvers/index.tsmodified+209−9View file
@@ -1,7 +1,12 @@
1-// The in-browser solvers shipped with the site. Each solver is a MATLAB
2-// function file (solver.m) implementing the interface documented in
1+// The solvers shipped with the site. Each solver is a MATLAB function
2+// file (solver.m) implementing the interface documented in
33 // docs/problems/laplace-dirichlet-2d.md. A submitted solver adds a
44 // directory here plus a manifest entry.
5+//
6+// Two manifest entries may share one solver.m: the "-mat" entries are the
7+// same file as their numbl twin, run in real MATLAB instead, which makes
8+// the pair of curves a measurement of the runtime rather than of the
9+// method. Such an entry sets sourceDir to the twin's directory.
510
611 export interface SolverManifest {
712 /** Short stable identifier used in results and URLs. */
@@ -18,6 +23,25 @@ export interface SolverManifest {
1823 runtime: "numbl" | "matlab";
1924 /** The resolution values a standard work-precision sweep runs. */
2025 sweepN: number[];
26+ /** Resolutions for instances that need a different range from sweepN,
27+ * keyed by instance id. A harder geometry can need several times the
28+ * resolution, and putting those values in every sweep would only make
29+ * the easy instances slow. */
30+ sweepNByInstance?: Record<string, number[]>;
31+ /** Directory under src/solvers/ holding this solver's solver.m, when it
32+ * is not the solver's own id: set by entries that share a file with
33+ * another entry. Both entries must then carry the same version. */
34+ sourceDir?: string;
35+}
36+
37+/** Directory under src/solvers/ holding a solver's solver.m. */
38+export function solverSourceDir(s: SolverManifest): string {
39+ return s.sourceDir ?? s.id;
40+}
41+
42+/** The resolutions a sweep of this solver runs on this instance. */
43+export function sweepNFor(s: SolverManifest, instanceId: string): number[] {
44+ return s.sweepNByInstance?.[instanceId] ?? s.sweepN;
2145 }
2246
2347 export const SOLVERS: SolverManifest[] = [
@@ -33,7 +57,29 @@ export const SOLVERS: SolverManifest[] = [
3357 "integral-equation methods. When the data's singularities sit inside " +
3458 "that curve, as on star-hard, convergence is lost: more charges keep " +
3559 "helping only until the system's ill-conditioning takes over, and the " +
36- "error settles near 1e-6 however far the sweep is pushed.",
60+ "error settles near 1e-6 however far the sweep is pushed. On " +
61+ "flower-15 the method fails outright, and for a reason worth stating: " +
62+ "a fixed offset of 0.3 from a boundary whose curvature radius is " +
63+ "0.014 is not a curve at all. The offset self-intersects and leaves " +
64+ "some 45 percent of the charges inside the domain, where they " +
65+ "destroy the harmonicity of the representation. A better MFS would " +
66+ "choose the offset from the local geometry rather than fixing it in " +
67+ "advance; this one does not, and the instance shows what that costs. " +
68+ "On square-corners it is instead the strongest method by a wide " +
69+ "margin, reaching 2e-13. Two things go its way there: the domain is " +
70+ "convex, so the offset charge curve is simple again, and the charges " +
71+ "stay 0.3 away from every evaluation point, so the near-corner " +
72+ "targets that defeat the integral-equation methods cost it nothing. " +
73+ "Its error at the target 0.005 inside a corner is the same 1e-13 as " +
74+ "in the bulk. On star-nearfield, whose targets come within 0.005 of " +
75+ "the boundary, it likewise pays nothing for them: charges 0.3 " +
76+ "outside the domain make the representation smooth all the way in, " +
77+ "so the error there is the conditioning-limited few times 1e-6 it " +
78+ "already reaches on star-hard (1.7e-6 in real MATLAB, 3.9e-6 in " +
79+ "numbl, the two differing because the error is set by rounding in a " +
80+ "badly conditioned solve). That makes it the cheapest route to 1e-6 " +
81+ "of the three solvers on that instance, and the one with nothing " +
82+ "left beyond it.",
3783 version: "1.0.0",
3884 backend: "cpu",
3985 runtime: "numbl",
@@ -47,11 +93,54 @@ export const SOLVERS: SolverManifest[] = [
4793 "density, discretized with the periodic trapezoid rule at n " +
4894 "boundary nodes. Converges geometrically for any smooth data, at a " +
4995 "rate set by how far the data continues analytically, so the hard " +
50- "instance costs more nodes rather than a lost method assumption.",
96+ "instance costs more nodes rather than a lost method assumption. On " +
97+ "flower-15 that price is about four times the nodes: a few hundred " +
98+ "to resolve fifteen lobes at all, and more still because the plain " +
99+ "quadrature is inaccurate for the evaluation points nearest the " +
100+ "boundary until the node spacing falls well below their distance to " +
101+ "it (0.0295 at the closest). Past that the geometric convergence " +
102+ "returns, and it is fast: on that instance the error falls from " +
103+ "5e-4 at n = 768 to 2e-7 at n = 1536. square-corners exposes the " +
104+ "method's real weakness instead, which is evaluation rather than " +
105+ "solution. At n = 768 the bulk error there is 7e-9, already " +
106+ "converged, while the error at the target 0.005 inside a corner is " +
107+ "4.6e-2: uncorrected quadrature for a target a fraction of a node " +
108+ "spacing from the boundary loses everything, and the reported error " +
109+ "follows exp(-2 pi delta / h) across the four target distances. " +
110+ "Reaching 1e-8 at those targets takes about 4000 nodes, where the " +
111+ "bulk alone would need 500. star-nearfield makes the same point " +
112+ "without the geometry, being the domain and data of star-hard with " +
113+ "32 targets 0.005 to 0.05 inside the boundary, and there the " +
114+ "exp(-2 pi delta / h) law is legible straight down the ladder: at " +
115+ "1024 nodes one and the same density is accurate to 7.5e-14 in the " +
116+ "bulk, 6e-14 at the target 0.05 in, 6.6e-6 at 0.02, 2.7e-3 at 0.01 " +
117+ "and 6.3e-2 at 0.005. The browser sweep, stopping at 1536 nodes, " +
118+ "reaches 1.2e-2; the MATLAB twin needs 4096 nodes and 0.86 s for " +
119+ "6.7e-6. We emphasize that this is a property of the uncorrected " +
120+ "quadrature rule and not of Nystrom methods: a near-field " +
121+ "correction, whether by kernel splitting or by a locally corrected " +
122+ "or expansion-based rule, removes it, and chunkie-dlp, which is the " +
123+ "same integral equation with one, reaches 1.8e-11 there.",
51124 version: "1.0.0",
52125 backend: "cpu",
53126 runtime: "numbl",
54127 sweepN: [16, 24, 32, 48, 64, 96, 128, 192, 256, 384, 512, 768],
128+ // flower-15 needs roughly four times the nodes: the boundary itself
129+ // costs about 200 to resolve, and the evaluation points sit close
130+ // enough to it that the plain quadrature stays inaccurate until the
131+ // node spacing is well under that distance.
132+ sweepNByInstance: {
133+ "flower-15": [32, 48, 64, 96, 128, 192, 256, 384, 512, 768, 1024, 1536],
134+ // The corner targets need a node spacing well below 0.005, which is
135+ // thousands of nodes. This list stops at 1536 because a numbl solve
136+ // at 2048 takes about five seconds, and a browser sweep has to stay
137+ // usable; nystrom-dlp-mat carries the tail.
138+ "square-corners": [64, 96, 128, 192, 256, 384, 512, 768, 1024, 1536],
139+ // Same story as square-corners, for the same reason: the closest
140+ // near-boundary target is 0.005 inside, and the coarse end of the
141+ // standard list is off the chart (relMax above 1 at n = 16).
142+ "star-nearfield": [64, 96, 128, 192, 256, 384, 512, 768, 1024, 1536],
143+ },
55144 },
56145 {
57146 id: "chunkie-dlp",
@@ -62,14 +151,125 @@ export const SOLVERS: SolverManifest[] = [
62151 "toolbox: n uniform 16th-order Gauss-Legendre panels, high-order " +
63152 "singular quadrature in the assembly, a direct dense solve, and " +
64153 "near-corrected evaluation of the potential. Runs in real MATLAB " +
65- "only: the command line invokes matlab -batch and fetches chunkie " +
66- "on first use, so its results appear here but cannot be rerun in " +
67- "the browser. chunkie's default quadrature tolerances cap the " +
68- "attainable accuracy near 1e-11.",
69- version: "2.0.0",
154+ "only: the command line invokes matlab -batch and installs chunkie " +
155+ "with the mip package manager on first use, so its results appear " +
156+ "here but cannot be rerun in the browser. chunkie's default " +
157+ "quadrature tolerances cap the attainable accuracy near 1e-11. Its " +
158+ "corrected quadrature buys a great deal of accuracy per node: on " +
159+ "flower-15 it is some 250 times more accurate than nystrom-dlp at " +
160+ "the same number of boundary points. Note, however, that this does " +
161+ "not make it the faster route to a given accuracy here. Each solve " +
162+ "costs roughly eight times as much per point (adaptive singular " +
163+ "quadrature in the assembly, near-corrections in the evaluation), " +
164+ "and against that the trapezoid rule's geometric convergence wins: " +
165+ "on flower-15 in the same MATLAB, nystrom-dlp-mat reaches 1e-9 in " +
166+ "0.2 s where chunkie needs 0.3 s, and goes on to 1e-14, which " +
167+ "chunkie's tolerances do not permit. The generality chunkie pays " +
168+ "for (corners, adaptive panels, fast algorithms) is barely " +
169+ "exercised by the star instances. The fast algorithms in particular " +
170+ "never get their chance. The solver runs with chunkie's defaults, " +
171+ "so the fmm2d library that mip installs alongside it accelerates " +
172+ "the evaluation, but with some 300 evaluation points there is " +
173+ "little there to accelerate: it is worth 1.2 times at the top of " +
174+ "the flower-15 sweep and nothing below. The two routes that would " +
175+ "replace the dense assembly and factorization, a chunkerflam fast " +
176+ "direct solve or GMRES on an FMM matvec, both break even only past " +
177+ "about 4000 boundary points, which on this problem is beyond where " +
178+ "the error has already stopped improving. square-corners exercises " +
179+ "part of it, and is one of the two instances where chunkie leads " +
180+ "the other integral-equation solver: its corrected evaluation " +
181+ "handles the " +
182+ "targets inside the corners, so what limits it is resolving a " +
183+ "corner of radius 0.014 with uniform panels rather than the " +
184+ "proximity of the targets. It reaches 2e-8 in 0.6 s where " +
185+ "nystrom-dlp-mat needs 0.85 s for 7e-9 and 0.23 s to manage only " +
186+ "1e-4. Note that both are beaten there by mfs, whose singularities " +
187+ "lie outside the domain and which therefore has no near-field problem " +
188+ "at all. The other instance, and the clearest case for the toolbox, " +
189+ "is star-nearfield. Its domain and data are star-hard's, so the " +
190+ "equation is identical and only the targets move, to within 0.005 " +
191+ "of the boundary; what is measured is the corrected near-field " +
192+ "evaluation by itself. chunkie loses nothing there: 1.8e-11 at 64 " +
193+ "panels in 0.11 s, against 2.3e-11 at 48 panels in 0.05 s on " +
194+ "star-hard, so the near targets cost it a third more panels and " +
195+ "about twice the time, and no accuracy. nystrom-dlp-mat, the same " +
196+ "integral equation with an uncorrected evaluation, gives 1.5e-10 on " +
197+ "star-hard in 0.017 s and 1.5e-1 here at those same 768 nodes, and " +
198+ "needs 4096 nodes and 0.86 s to bring the closest target to 6.7e-6: " +
199+ "eight times chunkie's time for five orders less accuracy. mfs-mat, " +
200+ "which as on square-corners has no near field to correct, is " +
201+ "stopped instead by its own conditioning at 1.7e-6, and is the " +
202+ "cheapest route to that accuracy (0.017 s against chunkie's " +
203+ "0.065 s for 2.4e-6). Past 1e-6, though, chunkie is the only one of " +
204+ "the three that goes anywhere at all.",
205+ version: "3.0.0",
70206 backend: "cpu",
71207 runtime: "matlab",
72208 sweepN: [2, 3, 4, 6, 8, 12, 16, 24, 32, 48],
209+ // Two panels cannot see a 15-lobe boundary at all, and the sweep has
210+ // to run out to 192 panels before the quadrature tolerance caps it.
211+ sweepNByInstance: {
212+ "flower-15": [8, 12, 16, 24, 32, 48, 64, 96, 128, 192],
213+ // Uniform panels have to get down to the corner radius of 0.014.
214+ "square-corners": [8, 12, 16, 24, 32, 48, 64, 96, 128, 192, 256],
215+ // Two panels past the standard list: on star-hard the quadrature
216+ // tolerance caps the error by n = 48, and the near-boundary targets
217+ // of star-nearfield take a little more resolution to reach the same
218+ // cap.
219+ "star-nearfield": [2, 3, 4, 6, 8, 12, 16, 24, 32, 48, 64, 96],
220+ },
221+ },
222+ {
223+ id: "mfs-mat",
224+ name: "Method of fundamental solutions (MATLAB)",
225+ description:
226+ "The mfs solver file, unchanged, executed by real MATLAB instead of " +
227+ "numbl. The method, the resolution sweep, and the source file are the " +
228+ "same, so the pair of curves measures the runtime rather than the " +
229+ "discretization. Where the method converges, the two agree to the " +
230+ "digits reported; in the conditioning-limited regime of star-hard " +
231+ "and star-nearfield they do not, since there the error is set by " +
232+ "rounding in a badly " +
233+ "conditioned solve and the two linear-algebra libraries round " +
234+ "differently. Runs from the command line only, which invokes " +
235+ "matlab -batch.",
236+ version: "1.0.0",
237+ backend: "cpu",
238+ runtime: "matlab",
239+ sweepN: [8, 12, 16, 24, 32, 48, 64, 96, 128, 192, 256, 384, 512, 768],
240+ sourceDir: "mfs",
241+ },
242+ {
243+ id: "nystrom-dlp-mat",
244+ name: "Nystrom double-layer BIE (MATLAB)",
245+ description:
246+ "The nystrom-dlp solver file, unchanged, executed by real MATLAB " +
247+ "instead of numbl, with the same resolution sweep except on " +
248+ "square-corners and star-nearfield, where being some fifteen times " +
249+ "faster lets it follow the near-field targets out to 4096 nodes, " +
250+ "past where a browser sweep would stay usable. Unlike " +
251+ "chunkie-dlp, which is a different discretization of the same " +
252+ "integral equation, nothing here changes the numerics: the systems " +
253+ "are well conditioned on every instance and the reported errors " +
254+ "reproduce those of nystrom-dlp, so what remains between the two " +
255+ "curves is the cost of the runtime. Runs from the command line only, " +
256+ "which invokes matlab -batch.",
257+ version: "1.0.0",
258+ backend: "cpu",
259+ runtime: "matlab",
260+ sweepN: [16, 24, 32, 48, 64, 96, 128, 192, 256, 384, 512, 768],
261+ sweepNByInstance: {
262+ "flower-15": [32, 48, 64, 96, 128, 192, 256, 384, 512, 768, 1024, 1536],
263+ // Real MATLAB is fast enough to follow the corner targets all the
264+ // way down, which the numbl twin's list stops short of.
265+ "square-corners": [
266+ 64, 96, 128, 192, 256, 384, 512, 768, 1024, 1536, 2048, 3072, 4096,
267+ ],
268+ "star-nearfield": [
269+ 64, 96, 128, 192, 256, 384, 512, 768, 1024, 1536, 2048, 3072, 4096,
270+ ],
271+ },
272+ sourceDir: "nystrom-dlp",
73273 },
74274 ];
75275
src/solvers/nystrom-dlp/solver.mmodified+7−4View file
@@ -9,10 +9,13 @@ function out = solver(prob, n)
99 % W is discretized with the periodic trapezoid rule at n equispaced
1010 % parameter nodes; the kernel is smooth on a smooth curve, with the
1111 % 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.
12+% at a rate set by how far the data g continues analytically. No
13+% correction is applied for targets close to the boundary, so the error
14+% of the evaluated potential at a target a distance delta inside behaves
15+% like exp(-2 pi delta / h) in the node spacing h (the close-evaluation
16+% problem). On most instances the evaluation points stay a modest
17+% distance inside and this costs nothing; on star-nearfield, whose
18+% closest target is 0.005 in, it is what limits the solver.
1619 %
1720 % n : number of boundary quadrature nodes.
1821
test/expected.tsadded+68−0View file
@@ -0,0 +1,68 @@
1+// What the shipped solvers must achieve, shared by both test suites
2+// (test/solver-test.ts through numbl, test/matlab-test.ts in real
3+// MATLAB). Keyed by the solver's source directory rather than its id, so
4+// that manifest entries sharing one solver.m — mfs and mfs-mat, and
5+// likewise nystrom-dlp — are held to the same expectation whichever
6+// runtime executes them.
7+
8+/** Best relMax the solver must reach somewhere in its full sweep. */
9+export const MUST_REACH: Record<string, Record<string, number>> = {
10+ // On flower-15 the MFS charge curve, a fixed 0.3 outside a boundary
11+ // whose curvature radius is 0.014, self-intersects and puts charges
12+ // inside the domain; the floor is only a sanity check that the solver
13+ // still returns something, and MUST_NOT_REACH below is the real test.
14+ mfs: {
15+ "disk-easy": 1e-12,
16+ "star-medium": 1e-12,
17+ "star-hard": 1e-4,
18+ "flower-15": 1e-2,
19+ // Convex geometry and charges far from every target: MFS is the best
20+ // method on the corner instance, so this floor is a real requirement.
21+ "square-corners": 1e-11,
22+ // Charges 0.3 out and targets 0.005 in: the charge curve overshoots
23+ // the sources, as on star-hard, and the near-boundary targets cost
24+ // nothing extra because the representation is smooth up to it.
25+ "star-nearfield": 1e-4,
26+ },
27+ "nystrom-dlp": {
28+ "disk-easy": 1e-10,
29+ "star-medium": 1e-10,
30+ "star-hard": 1e-8,
31+ "flower-15": 1e-6,
32+ // Loose because the numbl sweep stops at 1536 nodes, where the target
33+ // 0.005 inside a corner is still at 1e-3; the MATLAB twin, whose list
34+ // runs to 4096, gets to 7e-9 and passes the same floor.
35+ "square-corners": 1e-2,
36+ // Looser still: uncorrected quadrature at a target 0.005 from a
37+ // boundary resolved at h = 0.0059 (1536 nodes) is worth about 1e-2,
38+ // and the MATLAB twin at 4096 only reaches 7e-6.
39+ "star-nearfield": 5e-2,
40+ },
41+ // chunkie's default quadrature tolerances cap it near 1e-11, and its
42+ // sweep is checked on the extreme instances only.
43+ "chunkie-dlp": {
44+ "disk-easy": 1e-10,
45+ "star-hard": 1e-9,
46+ "flower-15": 1e-9,
47+ "square-corners": 1e-8,
48+ // The instance chunkie exists to win: its corrected near-field
49+ // quadrature has to hold the same accuracy at a target 0.005 from the
50+ // boundary that it reaches in the bulk.
51+ "star-nearfield": 1e-9,
52+ },
53+};
54+
55+/** Accuracy the solver must NOT reach. MFS is required to do badly on
56+ * both hard instances, for the two different reasons they exist: on
57+ * star-hard its charge curve lies beyond the data's singularities, and on
58+ * flower-15 the fixed offset self-intersects and puts charges inside the
59+ * domain. If either suddenly reached high accuracy, the instance would no
60+ * longer be testing what the spec says it tests. */
61+export const MUST_NOT_REACH: Record<string, Record<string, number>> = {
62+ mfs: { "star-hard": 1e-8, "flower-15": 1e-8 },
63+ // star-nearfield exists to measure near-field evaluation, and the plain
64+ // trapezoid rule has none: if this solver ever reached high accuracy at
65+ // a target 0.005 inside the boundary, either it acquired a near-field
66+ // correction or the instance stopped placing its targets there.
67+ "nystrom-dlp": { "star-nearfield": 1e-8 },
68+};
test/matlab-test.tsmodified+46−33View file
@@ -1,17 +1,24 @@
11 // Convergence test for the MATLAB-runtime solvers, run where real MATLAB
22 // exists (not in CI): npx tsx test/matlab-test.ts
33 // Exits quietly with a notice when no matlab is on the PATH.
4+//
5+// The expectations live in test/expected.ts and are keyed by source
6+// directory, so mfs-mat and nystrom-dlp-mat are held to exactly the same
7+// accuracy as their numbl twins in test/solver-test.ts: if a -mat curve
8+// misses it, the runtime and not the method is the suspect.
49
510 import { readFileSync } from "fs";
611 import { fileURLToPath } from "url";
712 import { dirname, join } from "path";
813 import { getInstance } from "../src/problems/laplace2d/spec";
9-import { getSolver } from "../src/solvers";
14+import { SOLVERS, solverSourceDir, sweepNFor } from "../src/solvers";
1015 import {
11- ensureChunkie,
1216 matlabAvailable,
17+ matlabSetup,
1318 runMatlabSweep,
1419 } from "../src/cli/matlabRun";
20+import { DEFAULT_TIMING } from "../src/harness/timing";
21+import { MUST_NOT_REACH, MUST_REACH } from "./expected";
1522
1623 if (!matlabAvailable()) {
1724 console.log("matlab not found on PATH; skipping MATLAB solver tests");
@@ -20,42 +27,48 @@ if (!matlabAvailable()) {
2027
2128 const root = join(dirname(fileURLToPath(import.meta.url)), "..");
2229 const read = (p: string) => readFileSync(join(root, p), "utf-8");
23-const sources = {
30+const base = {
2431 buildProblem: read("src/problems/laplace2d/matlab/build_problem.m"),
2532 bdata: read("src/problems/laplace2d/matlab/laplace2d_bdata.m"),
26- solver: read("src/solvers/chunkie-dlp/solver.m"),
27-};
28-
29-const mustReach: Record<string, number> = {
30- "disk-easy": 1e-10,
31- "star-hard": 1e-9,
3233 };
3334
3435 let failures = 0;
35-for (const [instId, reach] of Object.entries(mustReach)) {
36- console.log(`\n== ${instId} / chunkie-dlp (MATLAB)`);
37- console.log(" n relMax relL2 solve(s)");
38- let best = Infinity;
39- const { points, matlabVersion } = runMatlabSweep({
40- instance: getInstance(instId),
41- ns: getSolver("chunkie-dlp").sweepN,
42- repeats: 1,
43- sources,
44- setup: ensureChunkie(),
45- });
46- for (const p of points) {
47- best = Math.min(best, p.relMax);
48- console.log(
49- ` ${String(p.n).padStart(4)} ${p.relMax.toExponential(3)} ` +
50- `${p.relL2.toExponential(3)} ${p.solveSeconds.toFixed(4)}`
51- );
52- }
53- console.log(` (MATLAB ${matlabVersion})`);
54- if (best > reach) {
55- console.log(` FAIL: best relMax ${best.toExponential(2)} > ${reach}`);
56- failures++;
57- } else {
58- console.log(` ok (best relMax ${best.toExponential(2)})`);
36+for (const solver of SOLVERS.filter((s) => s.runtime === "matlab")) {
37+ const dir = solverSourceDir(solver);
38+ const sources = { ...base, solver: read(`src/solvers/${dir}/solver.m`) };
39+ const setup = matlabSetup(solver.id);
40+ for (const [instId, reach] of Object.entries(MUST_REACH[dir] ?? {})) {
41+ console.log(`\n== ${instId} / ${solver.id} (MATLAB)`);
42+ console.log(" n relMax relL2 solve(s)");
43+ let best = Infinity;
44+ const { points, matlabVersion } = runMatlabSweep({
45+ instance: getInstance(instId),
46+ ns: sweepNFor(solver, instId),
47+ timing: { ...DEFAULT_TIMING, minTimedRuns: 1, timeBudgetSeconds: 0 },
48+ sources,
49+ setup,
50+ });
51+ for (const p of points) {
52+ best = Math.min(best, p.relMax);
53+ console.log(
54+ ` ${String(p.n).padStart(4)} ${p.relMax.toExponential(3)} ` +
55+ `${p.relL2.toExponential(3)} ${p.solveSeconds.toFixed(4)}`
56+ );
57+ }
58+ console.log(` (MATLAB ${matlabVersion})`);
59+ const notReach = MUST_NOT_REACH[dir]?.[instId];
60+ if (best > reach) {
61+ console.log(` FAIL: best relMax ${best.toExponential(2)} > ${reach}`);
62+ failures++;
63+ } else if (notReach !== undefined && best < notReach) {
64+ console.log(
65+ ` FAIL: best relMax ${best.toExponential(2)} < ${notReach} ` +
66+ "(instance no longer defeats this method)"
67+ );
68+ failures++;
69+ } else {
70+ console.log(` ok (best relMax ${best.toExponential(2)})`);
71+ }
5972 }
6073 }
6174
test/solver-test.tsmodified+67−28View file
@@ -1,16 +1,26 @@
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
1+// Convergence test: run the numbl solvers on the official instances
2+// through numbl in node and check that errors behave as the theory says
3+// they should. The MATLAB-runtime solvers are covered by
4+// test/matlab-test.ts, run locally where MATLAB exists, against the same
5+// expectations in test/expected.ts.
6+// Run with: npx tsx test/solver-test.ts
47
5-import { readFileSync } from "fs";
8+import { existsSync, readFileSync } from "fs";
69 import { fileURLToPath } from "url";
710 import { dirname, join } from "path";
811 import { INSTANCES, getInstance } from "../src/problems/laplace2d/spec";
9-import { SOLVERS } from "../src/solvers";
12+import {
13+ SOLVERS,
14+ getSolver,
15+ solverSourceDir,
16+ type SolverManifest,
17+} from "../src/solvers";
1018 import { runPoint, type MatlabSources } from "../src/harness/runner";
1119 import { runSweep } from "../src/harness/sweep";
1220 import { setNumblFileIO } from "../src/harness/numblRun";
21+import { DEFAULT_TIMING } from "../src/harness/timing";
1322 import { NodeFileIOAdapter } from "../src/cli/nodeFileIO";
23+import { MUST_NOT_REACH, MUST_REACH } from "./expected";
1424
1525 // Solvers that mip-install packages (chunkie-dlp) need file I/O; in node
1626 // that is the curl-backed adapter.
@@ -23,27 +33,55 @@ const base = {
2333 buildProblem: read("src/problems/laplace2d/matlab/build_problem.m"),
2434 bdata: read("src/problems/laplace2d/matlab/laplace2d_bdata.m"),
2535 };
26-const solverSources: Record<string, MatlabSources> = {
27- mfs: { ...base, solver: read("src/solvers/mfs/solver.m") },
28- "nystrom-dlp": { ...base, solver: read("src/solvers/nystrom-dlp/solver.m") },
29-};
36+const solverSources = (solver: SolverManifest): MatlabSources => ({
37+ ...base,
38+ solver: read(`src/solvers/${solverSourceDir(solver)}/solver.m`),
39+});
3040
31-// Best relMax each solver must reach over its full sweep. On star-hard,
32-// MFS is additionally required NOT to do well: its charge curve lies
33-// beyond the data's singularities there, and if it suddenly reached high
34-// accuracy the instance would no longer be testing what the spec says.
35-const mustReach: Record<string, Record<string, number>> = {
36- mfs: { "disk-easy": 1e-12, "star-medium": 1e-12, "star-hard": 1e-4 },
37- "nystrom-dlp": { "disk-easy": 1e-10, "star-medium": 1e-10, "star-hard": 1e-8 },
38-};
39-const mustNotReach: Record<string, Record<string, number>> = {
40- mfs: { "star-hard": 1e-8 },
41-};
41+// The suite checks accuracy, not speed, so one timed run per point is
42+// enough; the committed results are what the full timing policy is for.
43+const TEST_TIMING = { ...DEFAULT_TIMING, minTimedRuns: 1, timeBudgetSeconds: 0 };
4244
4345 let failures = 0;
4446
45-// MATLAB-runtime solvers are covered by test/matlab-test.ts, run locally
46-// where MATLAB exists; this suite tests the numbl solvers.
47+// Registry consistency, checked for every entry including the ones this
48+// suite does not run: the solver file must exist, it must have stated
49+// expectations, and an entry that borrows another entry's solver.m must
50+// carry the same version, so that two results claiming the same solver
51+// version really did run the same code.
52+for (const s of SOLVERS) {
53+ const dir = solverSourceDir(s);
54+ const path = `src/solvers/${dir}/solver.m`;
55+ if (!existsSync(join(root, path))) {
56+ console.log(`FAIL: ${s.id} has no ${path}`);
57+ failures++;
58+ }
59+ const expected = MUST_REACH[dir];
60+ if (!expected) {
61+ console.log(`FAIL: ${s.id} has no entry in test/expected.ts`);
62+ failures++;
63+ } else if (s.runtime === "numbl") {
64+ for (const inst of INSTANCES) {
65+ if (expected[inst.id] === undefined) {
66+ console.log(`FAIL: ${s.id} has no expectation on ${inst.id}`);
67+ failures++;
68+ }
69+ }
70+ }
71+ const twin = s.sourceDir && SOLVERS.find((x) => x.id === s.sourceDir);
72+ if (s.sourceDir && !twin) {
73+ console.log(
74+ `FAIL: ${s.id} names sourceDir ${s.sourceDir}, which is not a solver`
75+ );
76+ failures++;
77+ } else if (twin && twin.version !== s.version) {
78+ console.log(
79+ `FAIL: ${s.id} v${s.version} shares solver.m with ${twin.id} v${twin.version}`
80+ );
81+ failures++;
82+ }
83+}
84+
4785 for (const inst of INSTANCES) {
4886 for (const solver of SOLVERS.filter((s) => s.runtime === "numbl")) {
4987 console.log(`\n== ${inst.id} / ${solver.id}`);
@@ -52,8 +90,8 @@ for (const inst of INSTANCES) {
5290 runSweep({
5391 instance: inst,
5492 solver,
55- sources: solverSources[solver.id],
56- repeats: 1,
93+ sources: solverSources(solver),
94+ timing: TEST_TIMING,
5795 onPoint: (p) => {
5896 best = Math.min(best, p.relMax);
5997 console.log(
@@ -62,8 +100,9 @@ for (const inst of INSTANCES) {
62100 );
63101 },
64102 });
65- const reach = mustReach[solver.id][inst.id];
66- const notReach = mustNotReach[solver.id]?.[inst.id];
103+ const dir = solverSourceDir(solver);
104+ const reach = MUST_REACH[dir][inst.id];
105+ const notReach = MUST_NOT_REACH[dir]?.[inst.id];
67106 if (best > reach) {
68107 console.log(` FAIL: best relMax ${best.toExponential(2)} > ${reach}`);
69108 failures++;
@@ -84,9 +123,9 @@ for (const inst of INSTANCES) {
84123 const p = runPoint({
85124 instance: getInstance("star-medium"),
86125 n: 64,
87- repeats: 1,
126+ timing: TEST_TIMING,
88127 wantGrid: true,
89- sources: solverSources["nystrom-dlp"],
128+ sources: solverSources(getSolver("nystrom-dlp")),
90129 });
91130 if (!p.uGrid || p.uGrid.length !== 200 * 200) {
92131 console.log(`\nFAIL: grid has ${p.uGrid?.length ?? 0} values, expected 40000`);