5127df5fastandaccurate: PDE solver benchmarks, laplace-dirichlet-2d with MFS and Nystrom DLP solversJeremy Magland 1# Problem: laplace-dirichlet-2d (version 1)
3This document is the canonical statement of the problem. Results refer to
4it by the pair (problem id, problem version); any change that could alter
5a measured number requires a version bump.
7## The problem
9Solve the Laplace equation
11 Δu = 0 in Ω, u = g on ∂Ω,
ad5dc23Three new instances, near-field target sets, and an adaptive timing policyJeremy Magland 13where Ω is the plane domain bounded by the curve
5127df5fastandaccurate: PDE solver benchmarks, laplace-dirichlet-2d with MFS and Nystrom DLP solversJeremy Magland 14
ad5dc23Three new instances, near-field target sets, and an adaptive timing policyJeremy Magland 15 x(t) = r(t) (cos t, sin t), t ∈ [0, 2π),
5127df5fastandaccurate: PDE solver benchmarks, laplace-dirichlet-2d with MFS and Nystrom DLP solversJeremy Magland 16
ad5dc23Three new instances, near-field target sets, and an adaptive timing policyJeremy Magland 17and g is Dirichlet data specified below. Two families of radius function
18are used, both analytic and both parameterized by the polar angle:
20 star: r(t) = 1 + a cos(k t)
21 rounded square: r(t) = (cos^p t + sin^p t)^(−1/p), p even
23The second traces the superellipse |x|^p + |y|^p = 1: four sides, straight
24to within about 1/p², meeting at four corners on the diagonals whose
25radius of curvature is about 1.4/p. Since p is even, cos^p t is a
26polynomial in cos t and the curve is analytic, so `curveDD` exists
27everywhere and a solver sees the same interface as for a star. A genuine
28corner, where the curve is only piecewise smooth and the layer density is
29singular, is a different problem and not this one; what this family gives
30is an arbitrarily *near* corner. The problem is posed in the
5127df5fastandaccurate: PDE solver benchmarks, laplace-dirichlet-2d with MFS and Nystrom DLP solversJeremy Magland 31continuum: a solver receives the curve and the data as functions, chooses
32its own discretization, and is scored on the values it returns at a fixed
33set of evaluation points.
35## The exact solution
59058e5Add chunkie-dlp solver: mip package support in the harness, curl-backed file I/O for nodeJeremy Magland 37The data g is manufactured from an exact harmonic function, a sum of
38three logarithmic point sources placed outside the domain:
5127df5fastandaccurate: PDE solver benchmarks, laplace-dirichlet-2d with MFS and Nystrom DLP solversJeremy Magland 39
40 u*(x) = Σ_j c_j log |x − s_j|, c = (1.0, −0.6, 0.8).
42Source s_j is the boundary point at parameter φ_j = 2π(j−1)/3 + 0.4
59058e5Add chunkie-dlp solver: mip package support in the harness, curl-backed file I/O for nodeJeremy Magland 43pushed a distance d along the outward unit normal, and g = u* restricted
44to ∂Ω. Since u* is harmonic in Ω, it is the unique solution, and errors
45are measured against it directly rather than against a reference
46computation. The distance d controls difficulty: u* continues
47harmonically only up to the sources, so the smaller d, the shorter the
48distance the data continues past the boundary, and methods whose
5127df5fastandaccurate: PDE solver benchmarks, laplace-dirichlet-2d with MFS and Nystrom DLP solversJeremy Magland 49representations assume a generous continuation lose it. A solver must
59058e5Add chunkie-dlp solver: mip package support in the harness, curl-backed file I/O for nodeJeremy Magland 50not use knowledge of the sources; they exist only to manufacture g.
5127df5fastandaccurate: PDE solver benchmarks, laplace-dirichlet-2d with MFS and Nystrom DLP solversJeremy Magland 51
52## Official instances
ad5dc23Three new instances, near-field target sets, and an adaptive timing policyJeremy Magland 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 |
63The hard instances are hard in different ways, which is the point of
64having several. On star-hard the geometry is mild and the data barely
65continues 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
67normal) and the boundary itself is the difficulty: fifteen lobes, with
68curvature up to 69, a curvature radius of 0.014 against a domain of size
691. Resolving it costs a few hundred boundary nodes before any method
70converges at all, and the fixed evaluation rule below puts the closest
71evaluation point 0.0295 from the boundary, so a method whose quadrature
72degrades near the boundary is held back further. Note that k = 15 is
73divisible by 3, so the three sources sit at equivalent phases of the lobe
74pattern; this is deliberate, and keeps the instance from favoring any
75particular angular sector.
77square-corners moves the difficulty from everywhere to somewhere. Its
78corners have radius 0.0142, one percent of the side, and its sides are
79straight to a thousandth: four short stretches of boundary demand all the
80resolution and the rest demands almost none. A method that can grade its
81discretization has something to exploit here that the star instances never
82offer, and the extra evaluation points described below sit inside the
83corners, at 0.005 to 0.05 from the boundary, so that near-field accuracy
84is measured where the geometry is worst rather than only where it is
85convenient.
87star-nearfield changes the question rather than the problem. Its domain
88and its data are exactly star-hard's, so the equation a solver has to
89solve is the same one; what differs is where the answer is wanted. Besides
90the standard points it carries 32 targets on the inward normal at eight
91parameters around the curve, at the same four distances the near-corner
92set uses, so the closest sits 0.005 from a boundary of perimeter 9.02. A
93method that evaluates its own representation with an uncorrected
94quadrature rule loses accuracy there for reasons that have nothing to do
95with how well it solved the equation: for the periodic trapezoid rule the
96error at a target a distance δ inside behaves like exp(−2πδ/h) in the node
97spacing h, so at n = 1024 nodes (h = 0.0088) one and the same computed
98density is accurate to 6e−14 at δ = 0.05, 6.6e−6 at 0.02, 2.7e−3 at 0.01
99and only 6.3e−2 at 0.005. The instance therefore measures near-field
100evaluation by itself, and the pair (star-hard, star-nearfield) isolates
101it: identical geometry, identical data, identical discretizations, one
102question changed. Note that a method whose representation is smooth up to
103the boundary, such as the MFS with its charges outside the domain, pays
104nothing here; the difficulty belongs to the representation and not to the
105problem.
5127df5fastandaccurate: PDE solver benchmarks, laplace-dirichlet-2d with MFS and Nystrom DLP solversJeremy Magland 106
107Committed results exist only at these instances, so every solver is
108compared on identical inputs. The parameters may be varied freely in the
109site's interactive views, but such runs are not recorded.
111## What is scored
ad5dc23Three new instances, near-field target sets, and an adaptive timing policyJeremy Magland 113Each instance fixes at least 289 evaluation points: 32 rays at angles
7fcfeb0Densify evaluation set to 289 points; evaluation-point overlay now a toggle, default offJeremy Magland 114θ_j = 2πj/32 + 0.13 (j = 0, …, 31), radial fractions
115ρ ∈ {0.1, 0.2, …, 0.9} along each ray (the point at fraction ρ on ray
5127df5fastandaccurate: PDE solver benchmarks, laplace-dirichlet-2d with MFS and Nystrom DLP solversJeremy Magland 116θ is ρ r(θ) (cos θ, sin θ)), ordered radius-major, plus the origin last.
118An instance with corners adds a **near-corner set** after those: for each
119corner in order (the diagonals θ_c = π/4 + jπ/2, j = 0, …, 3), the four
120points at distances δ ∈ {0.005, 0.01, 0.02, 0.05} inside the boundary
121along that diagonal, at radius r(θ_c) − δ. On square-corners this gives
12216 further points and 305 in all, and their true perpendicular distances
123to the boundary are 0.005, 0.010, 0.019 and 0.042, the smallest of them a
124third of the corner's own radius of curvature.
126An instance carrying the **near-boundary set** adds, after any near-corner
127set, four points per parameter: for each of the eight parameters
128t_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
130normal* at t_i, ordered with the parameter outer and the distance inner.
131On star-nearfield this gives 32 further points and 321 in all. Because
132every δ stays below the smallest radius of curvature
133there (0.0721), the distance from such a point to the curve is exactly δ,
134unlike the near-corner set whose radial offset makes the true distances a
135little smaller. The eight parameters are offset by 0.07 so that they land
136at no special phase of the lobe pattern, and 8 is coprime to k = 5, so
137they sample crests, flanks and valleys alike.
139The count is therefore per instance, and a solver must return one value
140per row of `prob.evalXY` rather than assume 289.
5127df5fastandaccurate: PDE solver benchmarks, laplace-dirichlet-2d with MFS and Nystrom DLP solversJeremy Magland 142The solver returns u at exactly these points. Reported errors are
144 relMax = max_i |u_i − u*_i| / max_i |u*_i|,
145 relL2 = ( Σ_i (u_i − u*_i)² / Σ_i u*_i² )^{1/2}.
ad5dc23Three new instances, near-field target sets, and an adaptive timing policyJeremy Magland 147relMax is the headline number; both are recorded. Note that the fraction
1480.9 bounds the *radial* distance, not the distance to the boundary. On a
149strongly wiggly instance the boundary runs nearly radially along the
150flanks of the lobes, so the perpendicular distance from an evaluation
151point to it is several times smaller than the radial one: 0.0295 at the
152closest point on flower-15, against 0.0464 on star-hard and 0.077 on
153star-medium. Close evaluation therefore enters on the wiggly instances,
154though not in its severe form. The severe form is star-nearfield, whose
155near-boundary set puts targets a fraction of a node spacing from the
156boundary. Adding it needed no version bump: it is a new instance, and no
157number measured on an existing one changes.
5127df5fastandaccurate: PDE solver benchmarks, laplace-dirichlet-2d with MFS and Nystrom DLP solversJeremy Magland 158
159## Timing protocol
161One run is a full call of the solver, including its own discretization,
2f05269Connect work-precision curves by resolution, not time; strengthen the timing protocolJeremy Magland 162assembly, solve, and evaluation. The harness performs two untimed warmup
ad5dc23Three new instances, near-field target sets, and an adaptive timing policyJeremy Magland 163runs (the first of which is recorded separately as the cold time), then
164timed runs, and reports the **fastest** of them as the solve time. The
165minimum is used rather than the mean or median because everything that
166interferes with a measurement (scheduling, other load, residual JIT
167compilation) only ever adds time, so the fastest run is the least
168contaminated estimate of the solver's own cost; every individual timing
169is recorded in the result file regardless.
171The number of timed runs is not fixed, because a fixed count times a
172cheap point badly: a solve of 0.1 ms is dominated by whatever else the
173machine was doing during those few hundred microseconds, and the minimum
174of five such samples still scatters by a factor of several from one sweep
175to the next, which made the low-resolution end of the curves noisy. So
176timed runs continue until they have accumulated a **time budget** (0.5 s
177by default), subject to a floor on the count (5) and a cap (50 runs).
178A solve of a fraction of a millisecond is therefore measured 50 times and
179an expensive one 5 times, and each point is sampled enough for its
180minimum to mean something. The cap also bounds how many timings a result
181file 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
184length of its `solveSecondsAll`. Results written before this policy
185record a fixed `protocol.timedRuns` instead. Where a
2f05269Connect work-precision curves by resolution, not time; strengthen the timing protocolJeremy Magland 186sweep shares one process across resolutions, as the MATLAB runner does,
187warmup runs precede the sweep as well, so that the first resolution does
188not absorb the session's one-time costs. All timing is tic/toc inside the
189solver's own runtime, so browser, node and MATLAB runs measure the same
190thing. Times from different machines are not comparable; every result
191records its environment, and comparisons across environments are the
192reader's responsibility.
194Note that solve time need not increase with the resolution parameter n. A
195solver whose cost is dominated by a per-target near-field computation can
196get *slower* as n decreases, since coarser panels put more targets in the
197near field. Work-precision curves are therefore drawn as parametric
198curves in n, and may double back in time; a point that is both slower and
199less accurate than another point on the same curve is simply dominated.
5127df5fastandaccurate: PDE solver benchmarks, laplace-dirichlet-2d with MFS and Nystrom DLP solversJeremy Magland 200
201## Solver interface
203A solver is a MATLAB function file
205 function out = solver(prob, n)
207where n is the solver's own resolution parameter (its meaning is the
208solver's choice; the standard sweep list is declared in the solver's
ad5dc23Three new instances, near-field target sets, and an adaptive timing policyJeremy Magland 209manifest, which may also declare a different list for an instance that
210needs a different range of resolutions) and prob is a struct with fields
5127df5fastandaccurate: PDE solver benchmarks, laplace-dirichlet-2d with MFS and Nystrom DLP solversJeremy Magland 211
212| field | meaning |
213|---|---|
214| `prob.curve` | `@(t) -> [x y]`, boundary point at parameter t (column vectors in, m×2 out) |
215| `prob.curveD` | first derivative of the curve with respect to t |
216| `prob.curveDD` | second derivative |
217| `prob.g` | `@(t) -> g`, Dirichlet data at boundary parameter t |
ad5dc23Three new instances, near-field target sets, and an adaptive timing policyJeremy Magland 218| `prob.evalXY` | m×2, the evaluation points (289; 305 with corners, 321 with the near-boundary set) |
5127df5fastandaccurate: PDE solver benchmarks, laplace-dirichlet-2d with MFS and Nystrom DLP solversJeremy Magland 219| `prob.vizXY` | m×2, visualization grid points (m = 0 when not requested) |
ad5dc23Three new instances, near-field target sets, and an adaptive timing policyJeremy Magland 221The return value is a struct: `out.uEval` (m×1, required) and
5127df5fastandaccurate: PDE solver benchmarks, laplace-dirichlet-2d with MFS and Nystrom DLP solversJeremy Magland 222`out.uGrid` (m×1; `[]` when `prob.vizXY` is empty). The solver must not
223reconstruct the sources analytically or otherwise special-case the known
224solution; submissions are reviewed for this.
acdea26Add mfs-gpu: the MFS on WebGPU, and a TypeScript form of the solver interfaceJeremy Magland 226### The TypeScript form
228A solver that cannot be a MATLAB file, one that runs on a GPU for
229instance, receives the same information as a plain object instead, built by
230`src/problems/laplace2d/problem.ts`:
232| field | meaning |
233|---|---|
234| `curve(t)` | `{x, y}`, the boundary point at parameter t |
235| `curveD(t)`, `curveDD(t)` | its first and second derivatives |
236| `g(t)` | Dirichlet data at boundary parameter t |
237| `evalXY` | `Float64Array`, `nEval` points as interleaved x, y |
238| `vizXY` | the visualization grid, interleaved, empty when not wanted |
240and returns `uEval`, and `uGrid` when the grid was asked for, as
241`Float64Array`. Everything else is the same in both forms: the same
242evaluation points in the same order, the same prohibition on
243reconstructing the sources, and the same timing protocol. Only the clock
244differs. A MATLAB solver is timed by tic/toc inside its own runtime; a GPU
245solver has no such clock, so it is timed on the host around a run that
246ends by awaiting the device, with the values read back before the clock
247stops. That is the same synchronization point MATLAB's tic/toc gives.
248Shader compilation and pipeline creation happen once per device rather
249than once per run, which puts them where numbl's JIT compilation already
250is: outside the timed runs, absorbed by the warmups.
252A note on what a GPU solver can and cannot compute here. WGSL has f32 and
253f16 and no double precision, and no extension in the WebGPU standard adds
254one. Emulating double in software (carrying a value as an unevaluated sum
255of two f32) is not a reliable way out either: it rests on error-free
256transformations such as `s = a + b; err = b - (s - a)`, which are exact
257only if the compiler evaluates them as written, and WGSL permits an
258implementation to use greater precision or to reassociate. So a WebGPU
259solver on this problem works in single precision, and for a method whose
260accuracy is limited by conditioning rather than by resolution that sets a
261ceiling the CPU does not have. `mfs-gpu` measures where that ceiling is.
5127df5fastandaccurate: PDE solver benchmarks, laplace-dirichlet-2d with MFS and Nystrom DLP solversJeremy Magland 263## Visualization grid
265When requested, `prob.vizXY` lists a 200×200 grid of points over the
ad5dc23Three new instances, near-field target sets, and an adaptive timing policyJeremy Magland 266bounding square [−R, R]², R = 1.05 max_t r(t), with flat index
5127df5fastandaccurate: PDE solver benchmarks, laplace-dirichlet-2d with MFS and Nystrom DLP solversJeremy Magland 267p = ix·200 + iy for x = xs[ix], y = xs[iy] (y varies fastest, MATLAB
268meshgrid column order). Points outside Ω are included and the viewer
269masks them; grid values are never scored.
271## Limitations
273The exact solution is smooth and free of boundary singularities, so this
274problem does not test corner handling, nonsmooth data, or interior
275sources (a nonzero right-hand side would exclude plain boundary-integral
59058e5Add chunkie-dlp solver: mip package support in the harness, curl-backed file I/O for nodeJeremy Magland 276methods; that belongs to a different problem). The domain family is
277star-shaped by construction, which some methods can exploit.
279Every boundary here is analytic, which bounds what the measured curves
280can show, and the instances differ in how much that bound bites. On the
281star instances the features are of one size everywhere: flower-15 is
282fifteen copies of the same lobe, so a method that refines adaptively has
283nothing to refine toward, and a global spectral rule, which must resolve
284the finest feature everywhere, loses nothing by doing so. There the
285geometric convergence of the global trapezoid rule is hard to beat.
286square-corners is the deliberate counterweight, with one localized feature
287and targets inside it, and the ranking there is different in kind rather
288than in degree. star-nearfield is a second counterweight, aimed at
289evaluation rather than at geometry: what orders the methods there is
290whether they correct their quadrature near the boundary at all, and a
291method that does not is held at 1e−2 where one that does reaches 1e−11.
293What remains untested is a genuine corner. The rounded square is smooth,
294so its layer density is smooth; a boundary with an actual vertex has a
295density that behaves like a power of the distance to it, needs graded
296refinement or a corner-specific quadrature, and can turn a well-posed
297second-kind equation into something that needs care. That is a different
298problem, not a harder instance of this one. Also untested is whether
299adaptive refinement pays: the panel solver here places its panels
300uniformly, so on square-corners it resolves the corner by refining the
301flat sides at the same time, and the benefit of grading is not measured.