5127df5fastandaccurate: PDE solver benchmarks, laplace-dirichlet-2d with MFS and Nystrom DLP solversJeremy Magland 1// Exact solution, geometry, and sampling rules for laplace-dirichlet-2d,
2// implemented independently of the MATLAB side (build_problem.m). The two
3// implementations of these simple formulas check each other: a solver can
4// only reach high accuracy if both agree.
ad5dc23Three new instances, near-field target sets, and an adaptive timing policyJeremy Magland 6import { hasCorners, hasNearBoundary, type Laplace2dInstance } from "./spec";
5127df5fastandaccurate: PDE solver benchmarks, laplace-dirichlet-2d with MFS and Nystrom DLP solversJeremy Magland 7
8export interface Source {
9 x: number;
10 y: number;
11 c: number;
12}
ad5dc23Three new instances, near-field target sets, and an adaptive timing policyJeremy Magland 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. */
5127df5fastandaccurate: PDE solver benchmarks, laplace-dirichlet-2d with MFS and Nystrom DLP solversJeremy Magland 18export function boundaryR(inst: Laplace2dInstance, t: number): number {
ad5dc23Three new instances, near-field target sets, and an adaptive timing policyJeremy Magland 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 }
5127df5fastandaccurate: PDE solver benchmarks, laplace-dirichlet-2d with MFS and Nystrom DLP solversJeremy Magland 23 return 1 + inst.a * Math.cos(inst.k * t);
24}
ad5dc23Three new instances, near-field target sets, and an adaptive timing policyJeremy Magland 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. */
28export 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}
40/** The largest radius the boundary reaches, which sets the view extent
41 * and the visualization grid. */
42export 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}
47/** Angles of the corners, for instances that have them: the four
48 * diagonals of the rounded square. */
49export 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}
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. */
59export const NEAR_DELTAS = [0.005, 0.01, 0.02, 0.05];
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. */
64export 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}
5127df5fastandaccurate: PDE solver benchmarks, laplace-dirichlet-2d with MFS and Nystrom DLP solversJeremy Magland 69/** Boundary point at parameter t. */
70export function boundaryPoint(inst: Laplace2dInstance, t: number) {
71 const r = boundaryR(inst, t);
72 return { x: r * Math.cos(t), y: r * Math.sin(t) };
73}
ad5dc23Three new instances, near-field target sets, and an adaptive timing policyJeremy Magland 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. */
79export 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}
59058e5Add chunkie-dlp solver: mip package support in the harness, curl-backed file I/O for nodeJeremy Magland 92/** The three exact-solution sources: boundary points at phi_j pushed a
5127df5fastandaccurate: PDE solver benchmarks, laplace-dirichlet-2d with MFS and Nystrom DLP solversJeremy Magland 93 * distance d along the outward normal, strengths [1.0, -0.6, 0.8]. */
94export function sources(inst: Laplace2dInstance): Source[] {
ad5dc23Three new instances, near-field target sets, and an adaptive timing policyJeremy Magland 95 const { d } = inst;
5127df5fastandaccurate: PDE solver benchmarks, laplace-dirichlet-2d with MFS and Nystrom DLP solversJeremy Magland 96 const strengths = [1.0, -0.6, 0.8];
97 return strengths.map((c, j) => {
59058e5Add chunkie-dlp solver: mip package support in the harness, curl-backed file I/O for nodeJeremy Magland 98 const phi = (2 * Math.PI * j) / 3 + 0.4;
ad5dc23Three new instances, near-field target sets, and an adaptive timing policyJeremy Magland 99 const r = boundaryR(inst, phi);
100 const dr = boundaryRD(inst, phi);
59058e5Add chunkie-dlp solver: mip package support in the harness, curl-backed file I/O for nodeJeremy Magland 101 const bx = r * Math.cos(phi);
102 const by = r * Math.sin(phi);
103 const dx = dr * Math.cos(phi) - r * Math.sin(phi);
104 const dy = dr * Math.sin(phi) + r * Math.cos(phi);
105 const sp = Math.hypot(dx, dy);
106 return { x: bx + (d * dy) / sp, y: by - (d * dx) / sp, c };
5127df5fastandaccurate: PDE solver benchmarks, laplace-dirichlet-2d with MFS and Nystrom DLP solversJeremy Magland 107 });
108}
59058e5Add chunkie-dlp solver: mip package support in the harness, curl-backed file I/O for nodeJeremy Magland 110/** Exact solution u(x, y) = sum_j c_j log|x - s_j|. */
5127df5fastandaccurate: PDE solver benchmarks, laplace-dirichlet-2d with MFS and Nystrom DLP solversJeremy Magland 111export function exactU(inst: Laplace2dInstance, x: number, y: number): number {
112 let u = 0;
113 for (const s of sources(inst)) {
114 u += s.c * 0.5 * Math.log((x - s.x) ** 2 + (y - s.y) ** 2);
115 }
116 return u;
117}
ad5dc23Three new instances, near-field target sets, and an adaptive timing policyJeremy Magland 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. */
5127df5fastandaccurate: PDE solver benchmarks, laplace-dirichlet-2d with MFS and Nystrom DLP solversJeremy Magland 127export function evalPoints(inst: Laplace2dInstance): { x: number; y: number }[] {
7fcfeb0Densify evaluation set to 289 points; evaluation-point overlay now a toggle, default offJeremy Magland 128 const rho = [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9];
5127df5fastandaccurate: PDE solver benchmarks, laplace-dirichlet-2d with MFS and Nystrom DLP solversJeremy Magland 129 const pts: { x: number; y: number }[] = [];
130 for (const r of rho) {
7fcfeb0Densify evaluation set to 289 points; evaluation-point overlay now a toggle, default offJeremy Magland 131 for (let j = 0; j < 32; j++) {
132 const th = (2 * Math.PI * j) / 32 + 0.13;
5127df5fastandaccurate: PDE solver benchmarks, laplace-dirichlet-2d with MFS and Nystrom DLP solversJeremy Magland 133 const rr = r * boundaryR(inst, th);
134 pts.push({ x: rr * Math.cos(th), y: rr * Math.sin(th) });
135 }
136 }
137 pts.push({ x: 0, y: 0 });
ad5dc23Three new instances, near-field target sets, and an adaptive timing policyJeremy Magland 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 }
5127df5fastandaccurate: PDE solver benchmarks, laplace-dirichlet-2d with MFS and Nystrom DLP solversJeremy Magland 149 return pts;
150}
152/** Exact solution at the evaluation points. */
153export function exactAtEvalPoints(inst: Laplace2dInstance): Float64Array {
154 const pts = evalPoints(inst);
155 const u = new Float64Array(pts.length);
156 pts.forEach((p, i) => (u[i] = exactU(inst, p.x, p.y)));
157 return u;
158}
160/** The visualization grid: ngrid x ngrid points over [-R, R]^2 with
ad5dc23Three new instances, near-field target sets, and an adaptive timing policyJeremy Magland 161 * R = 1.05 max_t r(t). Flat index p = ix * ngrid + iy with x = xs[ix],
5127df5fastandaccurate: PDE solver benchmarks, laplace-dirichlet-2d with MFS and Nystrom DLP solversJeremy Magland 162 * y = xs[iy] (y varies fastest), matching build_problem.m's meshgrid
163 * column order. */
164export const VIZ_NGRID = 200;
166export function vizGrid(inst: Laplace2dInstance) {
ad5dc23Three new instances, near-field target sets, and an adaptive timing policyJeremy Magland 167 const R = 1.05 * maxRadius(inst);
5127df5fastandaccurate: PDE solver benchmarks, laplace-dirichlet-2d with MFS and Nystrom DLP solversJeremy Magland 168 const xs = new Float64Array(VIZ_NGRID);
169 for (let i = 0; i < VIZ_NGRID; i++) {
170 xs[i] = -R + (2 * R * i) / (VIZ_NGRID - 1);
171 }
172 return { R, ngrid: VIZ_NGRID, xs };
173}
175/** Whether (x, y) is inside the domain (used only for display masking,
176 * never for scoring, so float tie-breaks at the boundary are harmless). */
177export function insideDomain(inst: Laplace2dInstance, x: number, y: number): boolean {
178 const rr = Math.hypot(x, y);
179 const th = Math.atan2(y, x);
180 return rr < boundaryR(inst, th);
181}
183/** Relative errors of numeric values against the exact solution at the
184 * evaluation points: max and L2, both relative to the exact values. */
185export function evalErrors(
186 inst: Laplace2dInstance,
187 uNum: ArrayLike<number>
188): { relMax: number; relL2: number } {
189 const uEx = exactAtEvalPoints(inst);
190 if (uNum.length !== uEx.length) {
191 throw new Error(`expected ${uEx.length} values, got ${uNum.length}`);
192 }
193 let maxDiff = 0;
194 let maxEx = 0;
195 let sumDiff2 = 0;
196 let sumEx2 = 0;
197 for (let i = 0; i < uEx.length; i++) {
198 const diff = Math.abs(uNum[i] - uEx[i]);
199 maxDiff = Math.max(maxDiff, diff);
200 maxEx = Math.max(maxEx, Math.abs(uEx[i]));
201 sumDiff2 += diff * diff;
202 sumEx2 += uEx[i] * uEx[i];
203 }
204 return {
205 relMax: maxDiff / maxEx,
206 relL2: Math.sqrt(sumDiff2 / sumEx2),
207 };
208}