088be4cAn interactive tutorial on the Dirichlet problem, accuracy against costJeremy Magland 1# CLAUDE.md
3Notes for anyone, human or agent, working in this repository. Read
4[README.md](README.md) first for what the tutorial contains.
6## The shape of the thing
8A Vite + React app that embeds [numbl](https://numbl.org) as a library
9(`numbl/browser`, `createNumblSession`). This is the seqlab / barycentric-rational
10pattern, not the `numbl-project.json` + site-viewer pattern: there is no project
11file and the deploy is a plain `vite build`.
13```
14src/
15 App.tsx sidebar, hash routing, chapter switching
16 ui/ Lab (one runnable section), ChapterView, Controls, Markdown
17 engine/ the numbl session wrapper
18 plot/ a small hand-rolled SVG layer, plus a canvas field plot
e5fc8c6An introduction and static figures, and the MATLAB written plainly againJeremy Magland 19 figures/ static diagrams, one .tsx per part, keyed in figures/index.ts
088be4cAn interactive tutorial on the Dirichlet problem, accuracy against costJeremy Magland 20 matlab/ the shared .m library, shown to the reader on its own page
21 chapters/
e5fc8c6An introduction and static figures, and the MATLAB written plainly againJeremy Magland 22 intro/ the introduction: front matter, no labs
088be4cAn interactive tutorial on the Dirichlet problem, accuracy against costJeremy Magland 23 onedim/driver.m the driver for chapters 1-4, plus their shared panels
24 twodim/driver.m the driver for chapters 5-7, plus their shared panels
25 chNN/ per-chapter text.md, editable .m script, chapter.ts
6292f12A benchmark chapter, and offline benchmarks to compare it againstJeremy Magland 26 ch11/panels.tsx also the offline-results upload for the benchmark plot
27bench/ offline ports of the two best methods for real MATLAB,
28 Python and Julia; the JSON they write (schema
29 "laplace-dirichlet-bench-1", documented in bench/README.md)
30 loads into chapter 11's work-precision plot
33A **chapter** is prose in one `text.md` with `<!--lab:id-->` markers where its labs
34go. A **lab** is a driver (fixed, hidden), an editable script (the method), a list of
35declarative controls, and a render function. Each lab owns one numbl session for as
36long as it is on screen: moving a control rewrites `params.json` and re-runs, while
37editing the script reboots the session, because the script is compiled into `main.m`
38at boot.
40`main.m` is `driver.m` + `"\n"` + the editable script, so the script's functions
41become **local functions of the driver** and the driver can call them. That is the
42whole of the contract; each driver's header comment states which functions its
43script must define.
e5fc8c6An introduction and static figures, and the MATLAB written plainly againJeremy Magland 45## Working with numbl
47**This tutorial needs numbl 0.4.19 or later.** It runs correctly on 0.4.18 and
48is drastically slower there: the L-shaped sweep in Chapter 9 takes about three
49minutes instead of six seconds, and Chapter 8 about six seconds instead of a
50fifth of one.
52The reason is worth recording. Most of what this section used to say was a list
53of ways to fall off numbl's JIT fast path, found by profiling this tutorial: a
54`switch` anywhere in a call chain, a struct parameter, `A(k,:) = row` with a
55loop index, a bare `tic;`, and the expensive one: the rule that nothing nested
56under a loop the JIT had declined was compiled either, so a single unsupported
57line could put a whole solver on the interpreter. They were fixed
58upstream in 0.4.19, and the drivers here were rewritten back into ordinary
59MATLAB once they were: the sweeps are plain `for` loops again. If you find a
60comment describing a recursion-instead-of-a-loop workaround, it is stale; the
61git history has the old shape.
63What is left is worth knowing.
65### Keep the *timed* loop type-stable
67numbl specialises a loop on the types flowing through it. A loop whose iterate
68is assigned from several different functions, or whose body indexes a cell
69array, is not something it can specialise, and it runs interpreted. That is
70usually harmless, since the expensive callee inside is still compiled, but when
71the loop being measured *is* the thing under test the overhead lands in the
72number on screen. Two consequences visible in `ch08/iterative.m`:
74- each iterative method has its own loop rather than one loop with a switch
75 inside it, and multigrid keeps its grid hierarchy in one flat vector rather
76 than a cell array;
77- residual and error histories accumulate as squares and are square-rooted
78 once at the end. A `sqrt` inside the loop costs about 60 per cent on the
79 sweep, because sqrt may return a complex value and the history array is then
80 no longer a plain real array.
82Neither is a correctness matter. Check with `%!numbl:assert_jit` inside a loop
83body (it throws if that unit ran interpreted) and `NUMBL_JIT_LOG=1`, which
84prints every JIT decline with its reason.
86### `jsondecode` of a JSON array of arrays gives a cell of column vectors
e5fc8c6An introduction and static figures, and the MATLAB written plainly againJeremy Magland 88MATLAB gives a matrix when the rows are the same length. The L-shaped driver
89unpacks `p.nsm` by hand.
91### `jsonencode` refuses complex values, and flattens single-row matrices
e5fc8c6An introduction and static figures, and the MATLAB written plainly againJeremy Magland 93`jsonencode([1 2 3])` is `[1,2,3]`, not `[[1,2,3]]`, so anything matrix-shaped
94goes through `lp_rows`, which converts it to a cell array of rows first. NaN and
95Inf cross as `null`, which is why the TypeScript side uses `Num = number | null`
96everywhere and every plot helper skips nulls. A value whose imaginary part is
97entirely zero counts as real and encodes fine.
99### Backslash densifies sparse matrices
e5fc8c6An introduction and static figures, and the MATLAB written plainly againJeremy Magland 101`A \ b` converts a sparse `A` to dense before calling LAPACK, so there is no
102sparse direct solver here. That is why Chapter 5 stops at $N = 48$ and why
103Chapters 6 to 10 are about methods that do not need one. The text says so where
104it matters.
106### Scalar loops are often faster than vectorised slices
e5fc8c6An introduction and static figures, and the MATLAB written plainly againJeremy Magland 108A hand-written double loop over a stencil beats the equivalent
109`A(2:end-1,...)` expression by about a factor of three, because the vectorised
110form allocates several temporaries per sweep. The kernels in
111`ch08/iterative.m` are written as loops for that reason, which also happens to
112make them readable.
1a8bb6aMake the cross-language benchmark measure the languagesJeremy Magland 114The same holds for building a dense kernel matrix, and there the measurement is
115worth keeping. Three shapes for the Nyström matrix, at n = 1896:
117| shape | numbl | MATLAB | Julia |
118|---|---|---|---|
119| scalar double loop | 50 ms | 9 ms | 7 ms |
120| implicit expansion, `x - x'` | 99 ms | 13 ms | 82 ms |
121| outer products, `x * ones(1,n)` | 117 ms | 62 ms | — |
123So the loop wins everywhere the tutorial runs, and `ch10/bie.m` and
124`ch11/pacman.m` both use it. The third row is what those files used to
125contain; the idiom is not required by numbl, which supports implicit expansion
126and is itself faster with it, and it was costing real MATLAB a factor of
127seven in that stage. If you find it anywhere else, it is worth an experiment.
128NumPy is the exception that proves the rule: it cannot run the loop at all, so
129the Python port in `bench/` assembles with whole-array expressions in place,
130and the chapter's prose makes that asymmetry part of the point.
088be4cAn interactive tutorial on the Dirichlet problem, accuracy against costJeremy Magland 132### The browser is slower than the CLI
134The CLI uses a native LAPACK addon; the browser falls back to a TypeScript
135implementation, and it announces which one on `console.error`. Dense solves are
e5fc8c6An introduction and static figures, and the MATLAB written plainly againJeremy Magland 136noticeably slower in the page than in `npm run test:matlab`. Nothing in the
137tutorial depends on the absolute numbers, but do not quote a CLI timing as if it
138were what a reader will see.
140## Measuring time honestly
142Timings come from `tic`/`toc` **inside** the MATLAB layer, around the solve only, so
143they exclude session boot, JSON encoding and rendering. Each method gets a warm-up
1a8bb6aMake the cross-language benchmark measure the languagesJeremy Magland 144call before its first timed solve, because the first call also pays for compilation.
145Below about 0.1 ms the measurement is dominated by fixed overhead and the
146work-precision curves flatten out on the left; that is the clock, not the method.
148Chapters 10 and 11, whose numbers are compared against other environments or quoted
149in the prose, keep the **fastest of several runs**: at least three whenever a solve
150is cheap enough to allow it, many more when it is very cheap, one when it is
151expensive, with the loop bounded to about 0.6 s per point.
153```matlab
154while rep == 0 || (rep < 50 && tot < 0.6 && (rep < 3 || tot < 0.1))
155```
157That rule is not cosmetic. A single measurement of a twenty-millisecond solve on
158this hardware varies by a factor of two with whatever else the machine is doing,
159which is enough to reorder neighbouring points on a work-precision plot and, in
160the offline benchmarks, to make one language look several times faster than
161another. The three scripts under `bench/` implement the same rule, and any new
162cross-environment comparison must too. The earlier chapters still use the older
163rule (repeat only when under 5 ms), which is fine because nothing outside those
164chapters is compared against their numbers.
166## Testing
168```
169npm run test:matlab # ~2 min: the .m layer through the numbl CLI
170npm run build && npm run test:browser [--shots]
171```
173Every assertion in `test:matlab` is a number stated in the prose: an observed order,
174an error at a given resolution, a ratio between two methods. If one fails, either the
175code or the text is wrong, and the text is not allowed to win by default.
177`test:browser` opens every chapter, presses run on the labs that wait to be asked,
178and checks that nothing throws and that the figures contain data. It cannot judge
179whether a picture is any good; `--shots` writes full-page screenshots into `tmp/` so
180that a human, or an agent with vision, can look.
e5fc8c6An introduction and static figures, and the MATLAB written plainly againJeremy Magland 182## Static figures
184`src/figures/` holds the diagrams that are drawn rather than computed: the problem
185statement, the stencils, the roundoff-floor V, the V-cycle, the re-entrant corner,
186the close-evaluation spike. They are addressed from the prose by `<!--fig:id-->`, the
187same way labs are addressed by `<!--lab:id-->`, and resolved through the registry in
188`figures/index.ts`. They exist because several chapters, Chapters 8 to 10 in
189particular, show nothing at all until the reader presses Run.
191They are plain inline SVG with no data behind them, written against a 716-unit-wide
192viewBox and `width="100%"`, using `plot/palette.ts` for colour so that a curve
193labelled `fd2` is the same blue everywhere. Two habits are worth keeping. Text is
194placed by hand, so a caption or a legend that fits at one size may collide at
195another: render and look. And the introduction's accuracy-against-cost picture is a
196schematic, which its caption says out loud; do not let it drift into looking like
197measured data.
199To look at them, point a headless browser at each chapter and screenshot
200`figure.static-figure`; `npm run test:browser --shots` writes whole pages, which is
201coarser but needs no extra script.
0b6ac9cBuild and publish the site from the vault's CI on every push to mainJeremy Magland 203## Deployment
e7743b1Rename .cofferdam to .mochiJeremy Magland 205`.mochi/workflows/deploy.yml` runs on every push to `main`: `npm ci`,
0b6ac9cBuild and publish the site from the vault's CI on every push to mainJeremy Magland 206`npm run build`, then `upload-pages-artifact` and `deploy-pages`, which the
207vault substitutes with its own implementations to publish the repository's
208site. Two constraints shape it. `actions/cache` is not implemented in the
209vault's CI, so `setup-node` must not be given `cache: npm`. And a vault site is
210served sandboxed, in an opaque origin, so `localStorage` and friends *throw*;
211the app uses none, and adding one would break the deployed page while working
212perfectly in `npm run dev`.
e7743b1Rename .cofferdam to .mochiJeremy Magland 214Watch a deploy with `mochi run list --limit 3 --json` and
215`mochi run view <n> --log`.
088be4cAn interactive tutorial on the Dirichlet problem, accuracy against costJeremy Magland 217## Plotting
219`src/plot/` is a small SVG layer with no chart library. `palette.ts` records which
220colour does which job and the validator run that justifies the set; if you add a
221series, re-run the validator rather than picking a hex that looks nice. Fields are
222drawn into a canvas at one pixel per grid cell and scaled, so a 512-by-512 grid costs
223no more to draw than a 16-by-16 one, and points outside the domain are passed as
224`null` and left as background.
226## Things deliberately not done
228- A sparse direct solver, for the reason above.
229- Graded meshes in Chapter 9. The chapter names local refinement and singular
230 enrichment as the standard remedies and then demonstrates a third, because the
231 third is the one that is short enough to read.
232- Close-evaluation quadrature in Chapter 10. The failure is shown and quantified; the
233 fix is a chapter of its own.
234- Neumann and mixed boundary conditions, anywhere.