1// The solvers shipped with the site. Each solver is a MATLAB function
2// file (solver.m) implementing the interface documented in
3// docs/problems/laplace-dirichlet-2d.md. A submitted solver adds a
4// directory here plus a manifest entry.
5//
6// 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.
11export interface SolverManifest {
12 /** Short stable identifier used in results and URLs. */
13 id: string;
14 name: string;
15 /** One-paragraph description of the method for the problem page. */
16 description: string;
17 /** Semantic version of the solver code; bump on any change that could
18 * alter results. */
19 version: string;
20 backend: "cpu" | "gpu";
21 /** What executes the solver: "numbl" solvers run in the browser and in
22 * the CLI; "matlab" solvers run only in real MATLAB via the CLI;
23 * "webgpu" solvers are TypeScript and WGSL rather than a MATLAB file and
24 * run wherever a WebGPU device can be had, which is the browser and, with
25 * the optional `webgpu` package, the CLI. */
26 runtime: "numbl" | "matlab" | "webgpu";
27 /** The resolution values a standard work-precision sweep runs. */
28 sweepN: number[];
29 /** Resolutions for instances that need a different range from sweepN,
30 * keyed by instance id. A harder geometry can need several times the
31 * resolution, and putting those values in every sweep would only make
32 * the easy instances slow. */
33 sweepNByInstance?: Record<string, number[]>;
34 /** Directory under src/solvers/ holding this solver's solver.m, when it
35 * is not the solver's own id: set by entries that share a file with
36 * another entry. Both entries must then carry the same version. */
37 sourceDir?: string;
38}
40/** Directory under src/solvers/ holding a solver's source. */
41export function solverSourceDir(s: SolverManifest): string {
42 return s.sourceDir ?? s.id;
43}
45/** The solver's source files, relative to src/solvers/. A MATLAB solver is
46 * one file; a WebGPU one is its TypeScript driver and the module that
47 * generates its WGSL. */
48export function solverFiles(s: SolverManifest): string[] {
49 const dir = solverSourceDir(s);
50 if (s.runtime === "webgpu") return [`${dir}/solver.ts`, `${dir}/wgsl.ts`];
51 return [`${dir}/solver.m`];
52}
54/** The resolutions a sweep of this solver runs on this instance. */
55export function sweepNFor(s: SolverManifest, instanceId: string): number[] {
56 return s.sweepNByInstance?.[instanceId] ?? s.sweepN;
57}
59export const SOLVERS: SolverManifest[] = [
60 {
61 id: "mfs",
62 name: "Method of fundamental solutions",
63 description:
64 "Represents the solution as n logarithmic point charges on a curve " +
65 "a fixed distance 0.3 outside the boundary, with strengths found by " +
66 "collocation at n boundary points. Converges geometrically when the " +
67 "data continues harmonically past the charge curve, reaching machine " +
68 "precision on the easier instances with far less work than the " +
69 "integral-equation methods. When the data's singularities sit inside " +
70 "that curve, as on star-hard, convergence is lost: more charges keep " +
71 "helping only until the system's ill-conditioning takes over, and the " +
72 "error settles near 1e-6 however far the sweep is pushed. On " +
73 "flower-15 the method fails outright, and for a reason worth stating: " +
74 "a fixed offset of 0.3 from a boundary whose curvature radius is " +
75 "0.014 is not a curve at all. The offset self-intersects and leaves " +
76 "some 45 percent of the charges inside the domain, where they " +
77 "destroy the harmonicity of the representation. A better MFS would " +
78 "choose the offset from the local geometry rather than fixing it in " +
79 "advance; this one does not, and the instance shows what that costs. " +
80 "On square-corners it is instead the strongest method by a wide " +
81 "margin, reaching 2e-13. Two things go its way there: the domain is " +
82 "convex, so the offset charge curve is simple again, and the charges " +
83 "stay 0.3 away from every evaluation point, so the near-corner " +
84 "targets that defeat the integral-equation methods cost it nothing. " +
85 "Its error at the target 0.005 inside a corner is the same 1e-13 as " +
86 "in the bulk. On star-nearfield, whose targets come within 0.005 of " +
87 "the boundary, it likewise pays nothing for them: charges 0.3 " +
88 "outside the domain make the representation smooth all the way in, " +
89 "so the error there is the conditioning-limited few times 1e-6 it " +
90 "already reaches on star-hard (1.7e-6 in real MATLAB, 3.9e-6 in " +
91 "numbl, the two differing because the error is set by rounding in a " +
92 "badly conditioned solve). That makes it the cheapest route to 1e-6 " +
93 "of the three solvers on that instance, and the one with nothing " +
94 "left beyond it.",
95 version: "1.0.0",
96 backend: "cpu",
97 runtime: "numbl",
98 sweepN: [8, 12, 16, 24, 32, 48, 64, 96, 128, 192, 256, 384, 512, 768],
99 },
100 {
101 id: "mfs-gpu",
102 name: "Method of fundamental solutions (WebGPU)",
103 description:
104 "The mfs method unchanged, with its assembly, its dense solve and " +
105 "its evaluation all on the GPU through WebGPU. It is the one solver " +
106 "here that is not a MATLAB file: TypeScript and WGSL, run on the " +
107 "page's own device in the browser and on Dawn from the command line. " +
108 "The solve is a right-looking LU with partial pivoting, three " +
109 "dispatches per column, with the right-hand side carried as an extra " +
110 "matrix column so that forward substitution disappears; at n = 768 " +
111 "that is some 3000 dispatches in a single submit. Everything is in " +
112 "f32, because WebGPU has no double precision and no extension in the " +
113 "standard adds one. For a method whose accuracy is set by " +
114 "conditioning rather than by resolution that is decisive: on " +
115 "disk-easy the error stops at 3.0e-7 where the same method in double " +
116 "reaches 9.5e-16, and likewise 4.0e-7 against 7.7e-16 on star-medium " +
117 "and 2.7e-6 against 2.3e-13 on square-corners. On the three " +
118 "instances where mfs is already conditioning-limited in double the " +
119 "gap narrows to two and a half or three orders (1.4e-3 against " +
120 "4.0e-6 on star-hard), since there both are losing to the same " +
121 "ill-conditioning and only the rate differs. Note where the ceiling " +
122 "arrives: n = 96 charges, which is far too small a dense solve to " +
123 "pay for a GPU. At that resolution this solver is 4.5 times slower " +
124 "than mfs in the same browser and ten times slower than mfs-mat in " +
125 "real MATLAB, and the dispatch overhead is plain at the bottom of " +
126 "the sweep, where 8 charges still cost a millisecond. It does become " +
127 "the faster of the two browser solvers past about n = 384, reaching " +
128 "6.6 times at n = 768, but f32 has taken its own error to 1e-5 by " +
129 "then. So the curve is dominated everywhere: at every accuracy this " +
130 "solver reaches, mfs on the CPU reaches it sooner. We emphasize that " +
131 "this is a statement about the method and the size of its systems " +
132 "rather than about the hardware. An O(n^3) solve does eventually pay " +
133 "for a GPU, and the same kernels would; the MFS simply cannot use an " +
134 "n that large, because single precision has ended its convergence " +
135 "long before. A method that stays accurate as n grows, which on this " +
136 "problem means the integral-equation solvers, is where a WebGPU " +
137 "backend would have something to win.",
138 version: "1.0.0",
139 backend: "gpu",
140 runtime: "webgpu",
141 // The same list as mfs, so the two curves land on the same resolutions
142 // and the pair isolates the backend. Most of the upper half is past
143 // where f32 stops improving, and is there to show both the ceiling and
144 // what the O(n^3) solve costs on either side.
145 sweepN: [8, 12, 16, 24, 32, 48, 64, 96, 128, 192, 256, 384, 512, 768],
146 },
147 {
148 id: "nystrom-dlp",
149 name: "Nystrom double-layer BIE",
150 description:
151 "Second-kind boundary integral equation for the double-layer " +
152 "density, discretized with the periodic trapezoid rule at n " +
153 "boundary nodes. Converges geometrically for any smooth data, at a " +
154 "rate set by how far the data continues analytically, so the hard " +
155 "instance costs more nodes rather than a lost method assumption. On " +
156 "flower-15 that price is about four times the nodes: a few hundred " +
157 "to resolve fifteen lobes at all, and more still because the plain " +
158 "quadrature is inaccurate for the evaluation points nearest the " +
159 "boundary until the node spacing falls well below their distance to " +
160 "it (0.0295 at the closest). Past that the geometric convergence " +
161 "returns, and it is fast: on that instance the error falls from " +
162 "5e-4 at n = 768 to 2e-7 at n = 1536. square-corners exposes the " +
163 "method's real weakness instead, which is evaluation rather than " +
164 "solution. At n = 768 the bulk error there is 7e-9, already " +
165 "converged, while the error at the target 0.005 inside a corner is " +
166 "4.6e-2: uncorrected quadrature for a target a fraction of a node " +
167 "spacing from the boundary loses everything, and the reported error " +
168 "follows exp(-2 pi delta / h) across the four target distances. " +
169 "Reaching 1e-8 at those targets takes about 4000 nodes, where the " +
170 "bulk alone would need 500. star-nearfield makes the same point " +
171 "without the geometry, being the domain and data of star-hard with " +
172 "32 targets 0.005 to 0.05 inside the boundary, and there the " +
173 "exp(-2 pi delta / h) law is legible straight down the ladder: at " +
174 "1024 nodes one and the same density is accurate to 7.5e-14 in the " +
175 "bulk, 6e-14 at the target 0.05 in, 6.6e-6 at 0.02, 2.7e-3 at 0.01 " +
176 "and 6.3e-2 at 0.005. The browser sweep, stopping at 1536 nodes, " +
177 "reaches 1.2e-2; the MATLAB twin needs 4096 nodes and 0.86 s for " +
178 "6.7e-6. We emphasize that this is a property of the uncorrected " +
179 "quadrature rule and not of Nystrom methods: a near-field " +
180 "correction, whether by kernel splitting or by a locally corrected " +
181 "or expansion-based rule, removes it, and chunkie-dlp, which is the " +
182 "same integral equation with one, reaches 1.8e-11 there.",
183 version: "1.0.0",
184 backend: "cpu",
185 runtime: "numbl",
186 sweepN: [16, 24, 32, 48, 64, 96, 128, 192, 256, 384, 512, 768],
187 // flower-15 needs roughly four times the nodes: the boundary itself
188 // costs about 200 to resolve, and the evaluation points sit close
189 // enough to it that the plain quadrature stays inaccurate until the
190 // node spacing is well under that distance.
191 sweepNByInstance: {
192 "flower-15": [32, 48, 64, 96, 128, 192, 256, 384, 512, 768, 1024, 1536],
193 // The corner targets need a node spacing well below 0.005, which is
194 // thousands of nodes. This list stops at 1536 because a numbl solve
195 // at 2048 takes about five seconds, and a browser sweep has to stay
196 // usable; nystrom-dlp-mat carries the tail.
197 "square-corners": [64, 96, 128, 192, 256, 384, 512, 768, 1024, 1536],
198 // Same story as square-corners, for the same reason: the closest
199 // near-boundary target is 0.005 inside, and the coarse end of the
200 // standard list is off the chart (relMax above 1 at n = 16).
201 "star-nearfield": [64, 96, 128, 192, 256, 384, 512, 768, 1024, 1536],
202 },
203 },
204 {
205 id: "chunkie-dlp",
206 name: "chunkie double-layer BIE",
207 description:
208 "The same second-kind double-layer formulation as nystrom-dlp, " +
209 "discretized by chunkie, a production MATLAB boundary-integral " +
210 "toolbox: n uniform 16th-order Gauss-Legendre panels, high-order " +
211 "singular quadrature in the assembly, a direct dense solve, and " +
212 "near-corrected evaluation of the potential. Runs in real MATLAB " +
213 "only: the command line invokes matlab -batch and installs chunkie " +
214 "with the mip package manager on first use, so its results appear " +
215 "here but cannot be rerun in the browser. chunkie's default " +
216 "quadrature tolerances cap the attainable accuracy near 1e-11. Its " +
217 "corrected quadrature buys a great deal of accuracy per node: on " +
218 "flower-15 it is some 250 times more accurate than nystrom-dlp at " +
219 "the same number of boundary points. Note, however, that this does " +
220 "not make it the faster route to a given accuracy here. Each solve " +
221 "costs roughly eight times as much per point (adaptive singular " +
222 "quadrature in the assembly, near-corrections in the evaluation), " +
223 "and against that the trapezoid rule's geometric convergence wins: " +
224 "on flower-15 in the same MATLAB, nystrom-dlp-mat reaches 1e-9 in " +
225 "0.2 s where chunkie needs 0.3 s, and goes on to 1e-14, which " +
226 "chunkie's tolerances do not permit. The generality chunkie pays " +
227 "for (corners, adaptive panels, fast algorithms) is barely " +
228 "exercised by the star instances. The fast algorithms in particular " +
229 "never get their chance. The solver runs with chunkie's defaults, " +
230 "so the fmm2d library that mip installs alongside it accelerates " +
231 "the evaluation, but with some 300 evaluation points there is " +
232 "little there to accelerate: it is worth 1.2 times at the top of " +
233 "the flower-15 sweep and nothing below. The two routes that would " +
234 "replace the dense assembly and factorization, a chunkerflam fast " +
235 "direct solve or GMRES on an FMM matvec, both break even only past " +
236 "about 4000 boundary points, which on this problem is beyond where " +
237 "the error has already stopped improving. square-corners exercises " +
238 "part of it, and is one of the two instances where chunkie leads " +
239 "the other integral-equation solver: its corrected evaluation " +
240 "handles the " +
241 "targets inside the corners, so what limits it is resolving a " +
242 "corner of radius 0.014 with uniform panels rather than the " +
243 "proximity of the targets. It reaches 2e-8 in 0.6 s where " +
244 "nystrom-dlp-mat needs 0.85 s for 7e-9 and 0.23 s to manage only " +
245 "1e-4. Note that both are beaten there by mfs, whose singularities " +
246 "lie outside the domain and which therefore has no near-field problem " +
247 "at all. The other instance, and the clearest case for the toolbox, " +
248 "is star-nearfield. Its domain and data are star-hard's, so the " +
249 "equation is identical and only the targets move, to within 0.005 " +
250 "of the boundary; what is measured is the corrected near-field " +
251 "evaluation by itself. chunkie loses nothing there: 1.8e-11 at 64 " +
252 "panels in 0.09 s, against 2.3e-11 at 48 panels in 0.05 s on " +
253 "star-hard, so the near targets cost it a third more panels and " +
254 "roughly twice the time, and no accuracy. nystrom-dlp-mat, the same " +
255 "integral equation with an uncorrected evaluation, gives 1.5e-10 on " +
256 "star-hard in 0.017 s and 1.5e-1 here at those same 768 nodes, and " +
257 "needs 4096 nodes and 0.86 s to bring the closest target to 6.7e-6: " +
258 "nine times chunkie's time for five orders less accuracy. mfs-mat, " +
259 "which as on square-corners has no near field to correct, is " +
260 "stopped instead by its own conditioning at 1.7e-6, and is the " +
261 "cheapest route to that accuracy (0.017 s against chunkie's " +
262 "0.042 s for 2.4e-6). Past 1e-6, though, chunkie is the only one of " +
263 "the three that goes anywhere at all.",
264 version: "3.0.0",
265 backend: "cpu",
266 runtime: "matlab",
267 sweepN: [2, 3, 4, 6, 8, 12, 16, 24, 32, 48],
268 // Two panels cannot see a 15-lobe boundary at all, and the sweep has
269 // to run out to 192 panels before the quadrature tolerance caps it.
270 sweepNByInstance: {
271 "flower-15": [8, 12, 16, 24, 32, 48, 64, 96, 128, 192],
272 // Uniform panels have to get down to the corner radius of 0.014.
273 "square-corners": [8, 12, 16, 24, 32, 48, 64, 96, 128, 192, 256],
274 // Two panels past the standard list: on star-hard the quadrature
275 // tolerance caps the error by n = 48, and the near-boundary targets
276 // of star-nearfield take a little more resolution to reach the same
277 // cap.
278 "star-nearfield": [2, 3, 4, 6, 8, 12, 16, 24, 32, 48, 64, 96],
279 },
280 },
281 {
282 id: "mfs-mat",
283 name: "Method of fundamental solutions (MATLAB)",
284 description:
285 "The mfs solver file, unchanged, executed by real MATLAB instead of " +
286 "numbl. The method, the resolution sweep, and the source file are the " +
287 "same, so the pair of curves measures the runtime rather than the " +
288 "discretization. Where the method converges, the two agree to the " +
289 "digits reported; in the conditioning-limited regime of star-hard " +
290 "and star-nearfield they do not, since there the error is set by " +
291 "rounding in a badly " +
292 "conditioned solve and the two linear-algebra libraries round " +
293 "differently. Runs from the command line only, which invokes " +
294 "matlab -batch.",
295 version: "1.0.0",
296 backend: "cpu",
297 runtime: "matlab",
298 sweepN: [8, 12, 16, 24, 32, 48, 64, 96, 128, 192, 256, 384, 512, 768],
299 sourceDir: "mfs",
300 },
301 {
302 id: "nystrom-dlp-mat",
303 name: "Nystrom double-layer BIE (MATLAB)",
304 description:
305 "The nystrom-dlp solver file, unchanged, executed by real MATLAB " +
306 "instead of numbl, with the same resolution sweep except on " +
307 "square-corners and star-nearfield, where being some fifteen times " +
308 "faster lets it follow the near-field targets out to 4096 nodes, " +
309 "past where a browser sweep would stay usable. Unlike " +
310 "chunkie-dlp, which is a different discretization of the same " +
311 "integral equation, nothing here changes the numerics: the systems " +
312 "are well conditioned on every instance and the reported errors " +
313 "reproduce those of nystrom-dlp, so what remains between the two " +
314 "curves is the cost of the runtime. Runs from the command line only, " +
315 "which invokes matlab -batch.",
316 version: "1.0.0",
317 backend: "cpu",
318 runtime: "matlab",
319 sweepN: [16, 24, 32, 48, 64, 96, 128, 192, 256, 384, 512, 768],
320 sweepNByInstance: {
321 "flower-15": [32, 48, 64, 96, 128, 192, 256, 384, 512, 768, 1024, 1536],
322 // Real MATLAB is fast enough to follow the corner targets all the
323 // way down, which the numbl twin's list stops short of.
324 "square-corners": [
325 64, 96, 128, 192, 256, 384, 512, 768, 1024, 1536, 2048, 3072, 4096,
326 ],
327 "star-nearfield": [
328 64, 96, 128, 192, 256, 384, 512, 768, 1024, 1536, 2048, 3072, 4096,
329 ],
330 },
331 sourceDir: "nystrom-dlp",
332 },
333];
335export function getSolver(id: string): SolverManifest {
336 const s = SOLVERS.find((x) => x.id === id);
337 if (!s) throw new Error(`Unknown solver: ${id}`);
338 return s;
339}