1% Driver for the benchmark chapter: one hard problem, every method.
2%
3% The editable script supplies
4%
5% names = method_names()
6% v = solve_pac(mid, n, g, px, py)
7% d = dof_pac(mid, n)
8%
9% and must return the computed solution at the points (px, py), whatever grid,
10% basis, or boundary discretisation it used internally. n is the method's own
11% resolution parameter: intervals per unit length for the difference scheme,
12% basis terms for the expansion, refinement depth for the panels, points per
13% direction for the collocation. dof_pac reports the number of unknowns behind
14% that parameter, which the page shows so that the parameters stay comparable.
15%
16% g(x, y) is the Dirichlet data: the exact solution restricted to the
17% boundary. The domain is the Pac-Man
18%
19% 0 < r < 1, 0 < theta < 3*pi/2,
20%
21% and the exact solutions are those of the L-shaped chapter (lp_lexact), whose
22% corner sits at the same origin with the same interior angle.
23p = lp_params();
24kind = p.kind;
25names = method_names();
26nm = numel(names);
27g = @(xx, yy) lp_lexact(kind, xx, yy);
29out = struct();
30out.names = names;
32% the boundary, drawn over the field pictures: the edge theta = 0, the arc,
33% and the edge theta = 3*pi/2, walked as one closed polyline
34ta = (0:120)' * (3*pi/2) / 120;
35e1 = linspace(0, 1, 15)';
36e2 = linspace(-1, 0, 15)';
37out.cx = [e1; cos(ta); zeros(15, 1)]';
38out.cy = [zeros(15, 1); sin(ta); e2]';
40if strcmp(p.mode, 'field')
41 % a picture on a fixed grid of spacing 1/32, masked to the open domain;
42 % pixels on the two straight edges are left out so that the boundary
43 % methods are never evaluated exactly on the curve they integrate over
44 m = 34;
45 xs = (-m:m)' / 32;
46 [X, Y] = meshgrid(xs, xs);
47 R2 = X.^2 + Y.^2;
48 onedge = (abs(Y) < 1e-12 & X > -1e-12) | (abs(X) < 1e-12 & Y < 1e-12);
49 inside = R2 < 1 - 1e-9 & ~(X > 1e-12 & Y < -1e-12) & ~onedge;
50 px = X(inside);
51 py = Y(inside);
52 t0 = tic;
53 v = solve_pac(p.mid, p.n, g, px, py);
54 out.ms = toc(t0) * 1000;
55 ve = g(px, py);
56 out.maxerr = max(abs(v - ve));
57 out.unknowns = dof_pac(p.mid, p.n);
59 U = NaN(2*m+1, 2*m+1);
60 E = NaN(2*m+1, 2*m+1);
61 U(inside) = v;
62 E(inside) = abs(v - ve);
63 out.x = xs';
64 out.y = xs';
65 out.U = lp_rows(U);
66 out.Err = lp_rows(E);
67else
68 % errors at the nodes of a fixed lattice of spacing 1/8, kept at least
69 % 0.15 from the arc and 0.25 from the straight edges, and reaching within
70 % 0.125 of the re-entrant corner (but excluding the corner itself, which
71 % is a boundary point)
72 [X, Y] = meshgrid((-8:8)' / 8);
73 R2 = X.^2 + Y.^2;
74 keep = R2 <= 0.7225 + 1e-9 & R2 > 1e-12 & ~(X > 1e-12 & Y < -1e-12) ...
75 & ~(X > 1e-12 & Y < 0.25 - 1e-9) & ~(Y < -1e-12 & X > -0.25 + 1e-9);
76 px = X(keep);
77 py = Y(keep);
78 ve = g(px, py);
80 % Each method gets its own list of resolutions. numbl's jsondecode turns
81 % a JSON array of arrays into a cell of column vectors where MATLAB gives
82 % a matrix; the offline harness runs this same file in real MATLAB, so
83 % both shapes are accepted.
84 NSc = p.nsm;
85 if iscell(NSc)
86 nn = numel(NSc{1});
87 NS = zeros(nm, nn);
88 for k = 1:nm
89 row = NSc{k};
90 NS(k, :) = row(:)';
91 end
92 else
93 NS = NSc;
94 nn = size(NS, 2);
95 end
96 E = NaN(nm, nn);
97 T = NaN(nm, nn);
98 Q = NaN(nm, nn);
99 DOF = NaN(nm, nn);
100 for k = 1:nm
101 % a ladder of zeros skips the method entirely, which is how the
102 % offline harness runs a subset of the methods through this driver
103 nw = 0;
104 for j = 1:nn
105 if isfinite(NS(k, j)) && NS(k, j) > 0
106 nw = NS(k, j);
107 break;
108 end
109 end
110 if nw == 0
111 continue;
112 end
113 vw = solve_pac(k, nw, g, px(1), py(1)); % warm-up
114 for j = 1:nn
115 n = NS(k, j);
116 if ~isfinite(n) || n <= 0
117 continue;
118 end
119 % Keep the fastest of several runs. Three runs at least whenever a solve
120 % is cheap enough to allow it, many more when it is very cheap, and one
121 % when it is expensive: the budget below caps the wall clock per point.
122 % The minimum over runs is what makes these numbers comparable at all,
123 % here and in the offline scripts, which apply the same rule. A single
124 % measurement of a fifty-millisecond solve reports whatever else the
125 % machine happened to be doing at the time.
126 el = Inf;
127 tot = 0;
128 rep = 0;
129 while rep == 0 || (rep < 50 && tot < 0.6 && (rep < 3 || tot < 0.1))
130 t0 = tic;
131 v = solve_pac(k, n, g, px, py);
132 dt = toc(t0);
133 el = min(el, dt);
134 tot = tot + dt;
135 rep = rep + 1;
136 end
137 T(k, j) = el * 1000;
138 E(k, j) = max(abs(v - ve));
139 DOF(k, j) = dof_pac(k, n);
140 end
141 Q(k, :) = lp_orders(E(k, :), NS(k, :));
142 end
143 out.ns = NS(1, :);
144 out.NS = lp_rows(NS);
145 out.DOF = lp_rows(DOF);
146 out.E = lp_rows(E);
147 out.T = lp_rows(T);
148 out.orders = lp_rows(Q);
149end
151lp_write(out);