Elliptical plates, a WebGPU time loop for the explicit method, and WAV export
Plate shape: 'rectangle' or 'ellipse' (a disk when Lx = Ly) within the same
Lx by Ly box. plate_domain.m marks the plate's footprint on the grid as a
staircase straddling the true edge, and plate_matrices.m is now a general
assembly on that footprint that reproduces the previous rectangle operators
exactly. Every edge condition works on the ellipse in both finite-difference
methods; the modal method refuses it. A simply supported curved edge gets the
curvature term the staircase loses (the plate paradox), so the disk converges
to the Bessel-function frequencies; tests check clamped, simply supported and
free disks and a 2:1 ellipse. Fields are NaN off the plate and drawn clipped
to the outline. Presets: a 300 mm disk, a 300 x 200 ellipse, the cymbal
stand-in as a disk.
Explicit finite differences: plate_fd.m split into setup, loop and output,
with a WebGPU port of the loop (gpu/fdLoop.ts) that runs in the worker when
the option is on, in single precision, with the setup and output still in
numbl. The sparse form of the loop (free edges, any ellipse) uses a per-node
keep mask. A Node test compares the port with the numbl loop through Dawn.
Audio: the radiated sound can be saved as a WAV file.
58 changed files+3213−440
README.mdmodified+19−4View file
@@ -1,6 +1,6 @@
11 # Plate Mechanics Lab
22
3-An interactive lab for the mechanics of thin metal plates. You specify geometry, material, edge conditions and an excitation; the app solves the damped Kirchhoff-Love equation and shows displacement, surface stress (von Mises and principal) over time, a probe time series and spectrum, and the radiated sound at a listener point, which you can play.
3+An interactive lab for the mechanics of thin metal plates. You specify geometry (a rectangle, or an ellipse such as a disk), material, edge conditions and an excitation; the app solves the damped Kirchhoff-Love equation and shows displacement, surface stress (von Mises and principal) over time, a probe time series and spectrum, and the radiated sound at a listener point, which you can play or save as a WAV file.
44
55 The numerics run in [numbl](https://github.com/flatironinstitute/numbl), a MATLAB-syntax environment executing in the browser. The MATLAB source shown on screen is the source that ran: edit any file and rerun. No physics is computed in JavaScript; the host only draws arrays and plays samples.
66
@@ -21,11 +21,23 @@ The app and the tutorials are peers over the same core and components.
2121 Three methods solve the same problem specification, and the app can run any two and compare them directly:
2222
2323 - **Modal superposition** (`plate_modal.m`): closed-form modes of the simply supported rectangle, each integrated exactly in time for piecewise-linear forcing. No time-step limit; the approximation is the modal truncation.
24-- **Explicit finite differences** (`plate_fd.m`): second-order scheme with the 13-point biharmonic stencil; simply supported, clamped or free edges (free edges use the energy-method operators of `plate_matrices.m`). Stable only below the CFL bound from `plate_fd_plan.m`; the slider goes past it on purpose.
24+- **Explicit finite differences** (`plate_fd.m`): second-order scheme with the 13-point biharmonic stencil; simply supported, clamped or free edges, or free edges pinned at the centre like a cymbal on its stand (the last two use the energy-method operators of `plate_matrices.m`). Stable only below the CFL bound from `plate_fd_plan.m`; the slider goes past it on purpose. The work is split into `plate_fd_setup.m`, `plate_fd_loop.m` and `plate_fd_output.m`, and the loop has a WebGPU port (below).
2525 - **Implicit finite differences** (`plate_implicit.m`): the same stencils assembled as sparse matrices (`plate_matrices.m`), stiffness averaged over three time levels (theta scheme, theta = 1/4 is the trapezoidal rule). One sparse Cholesky factorization, two triangular solves per step, no stability limit; the step is chosen for accuracy. Needs numbl 0.4.21 or later for its sparse direct solvers.
2626
2727 A method is a data record (`packages/core/src/methods/*.ts`) plus a `.m` file. Adding one means adding both; nothing else in core changes.
2828
29+## Shapes
30+
31+The plate is a rectangle or the ellipse inscribed in the same Lx by Ly box (a disk when Lx = Ly). The finite-difference methods embed the ellipse in the rectangular grid: `plate_domain.m` marks the nodes within half a spacing of the inside, so the edge is a staircase that straddles the true curve, and `plate_matrices.m` assembles every edge condition on that footprint as a sparse operator (for the rectangle it reproduces the old operators exactly). Supported and clamped edges hold the outermost ring of nodes at zero, the clamped ring with the mirrored Laplacian that gives the 7/dx^4 coefficient on a straight edge; free edges use the energy method on every node of the plate. The geometry, the edge conditions and the stress along the edge are then first order in the spacing, against second order inside, and the run notes say so.
32+
33+A simply supported curved edge needs one more term. With w = 0 on the edge the bending energy is D/2 [∫(∇²w)² dA − (1 − ν)∮ κ w_n² ds], with κ the curvature of the edge; a polygon or a staircase has κ = 0, so a scheme that only imposes w = 0 converges to a plate with zero Laplacian on the edge rather than zero moment (the Babuška, or Sapondzhyan, paradox). For a disk that puts the fundamental at λ² = 5.78 instead of 4.94 (ν = 0.3). `plate_matrices.m` adds the boundary term with the curvature of the true ellipse at each edge node, and the test suite checks the disk against the Bessel-function frequency equations: clamped 10.216, simply supported 4.935, free 5.358 and 9.003, all within a few percent on a 25-node grid and converging under refinement. The modal method has closed-form modes only for the simply supported rectangle and refuses the ellipse.
34+
35+## The GPU option
36+
37+The explicit scheme is many cheap steps, which is what a GPU does well, so its time loop has a WebGPU port: `packages/core/src/gpu/fdLoop.ts` is `plate_fd_loop.m` as compute shaders, three dispatches per step for a rectangle with supported or clamped edges (Laplacian, update, reductions for the probe, the audio bins and the blow-up test), three with the sparse operators of an ellipse, and five when the rigid-body projection of a free plate is added. A checkbox in the explicit method's settings turns it on when the browser exposes WebGPU. The setup (`plate_fd_setup.m`) and the output (`plate_fd_output.m`, stresses, sound, spectra) still run in numbl, unchanged and in double precision; the worker hands the setup struct to the GPU and the loop's record back. The generated script says which way it ran, the run details name the adapter, and editing `plate_fd_loop.m` sends the loop back to numbl, where the edit is live. This follows the pattern of the cymbal2 and turing-surface projects: the `.m` files remain the source of truth, and the port must compute the same record from the same struct.
38+
39+The GPU works in single precision. The state is carried as (w, w - w_prev) rather than (w, w_prev), the same arithmetic rearranged so that the rounding of w does not feed back into the velocity; measured against the double-precision loop over 50 ms on the 61 x 41 grid, the displacement differs by a relative L2 error of about 1e-6 and the stress and the sound by about 1e-4. The cost is about 20 µs per step on an integrated Intel GPU, nearly independent of the grid at these sizes: the 61 x 41 grid for 50 ms (27,000 steps) takes 0.8 s in total against 22 s in numbl, and a 121 x 81 grid (107,000 steps) takes 3 s. The implicit method is not ported: its step is two sparse triangular solves, which have no parallelism to offer at this size.
40+
2941 ## The live app
3042
3143 The lab is deployed at <https://plate-mechanics-sim.vault1-sites.magland.org/>. `.mochi/workflows/deploy.yml`
@@ -44,6 +56,8 @@ npm run typecheck
4456 npm run build
4557 ```
4658
59+`packages/core/test/gpu.test.ts` runs the WebGPU loop from Node against desktop WebGPU through the optional `webgpu` package (prebuilt Google Dawn) and compares it with the numbl loop for every edge condition on both shapes; `npm i --no-save webgpu` installs it. Without the package or a GPU adapter the suite skips itself.
60+
4761 The `.m` files also run unchanged from the numbl CLI, Octave or MATLAB:
4862
4963 ```bash
@@ -53,8 +67,9 @@ npx numbl run --path packages/core/src/m my_script.m
5367 ## Honest limits
5468
5569 - Kirchhoff-Love thin-plate theory: no shear deformation or rotary inertia, so thick plates and very high modes are wrong.
56-- Edge conditions are idealizations. Free edges are handled by a discrete bending energy whose natural boundary conditions are first-order accurate at the edge, and the rigid-body motion of a free plate is removed rather than modelled, as the supports of a real free-free rig would do.
70+- Edge conditions are idealizations. Free edges are handled by a discrete bending energy whose natural boundary conditions are first-order accurate at the edge, and the rigid-body motion of a free plate is removed rather than modelled, as the supports of a real free-free rig would do. The centre-pinned plate is held at one grid node, so the stress right at the support is resolution-dependent and the rocking a real stand's washers would restrain is removed instead.
71+- An ellipse is a staircase on the grid. Its edge is first order in the grid spacing, the stress on the outermost ring of nodes is copied from the nearest node with a full stencil (so the peak stress of a clamped edge sits about one spacing inside), and degenerate mode pairs of a disk are split slightly by the grid's anisotropy. Refine the grid to check.
5772 - Damping is a two-parameter model, chosen, not measured.
5873 - Stress is a second derivative of the displacement and converges more slowly than displacement in both methods; near the strike it depends on the load footprint.
5974 - Radiated sound is Rayleigh's integral in an infinite baffle: an estimate of the waveform, not a calibrated level.
60-- Explicit finite differences cost about 0.3 µs per cell-step in the browser, and the stable step scales with dx². The default 41 x 27 grid gives 20 ms in a few seconds; refining to 61 x 41 costs about five times more. The implicit scheme removes the step limit (about 1 ms per step on 61 x 41) at the price of detuning the highest modes; long simulated times still belong to the modal method.
75+- Explicit finite differences cost about 0.3 µs per cell-step in numbl, and the stable step scales with dx². The default 41 x 27 grid gives 20 ms in a few seconds; refining to 61 x 41 costs about five times more. The GPU option brings this to about 20 µs per step whatever the grid, in single precision. The implicit scheme removes the step limit (about 1 ms per step on 61 x 41) at the price of detuning the highest modes; long simulated times still belong to the modal method.
apps/lab/src/App.tsxmodified+49−22View file
@@ -4,6 +4,7 @@ import {
44 LIBRARY_FILES,
55 METHODS,
66 buildCompareBundle,
7+ buildGpuStages,
78 buildPlanBundle,
89 buildRunBundle,
910 defaultSpec,
@@ -29,6 +30,7 @@ import {
2930 RunStatus,
3031 useRunner,
3132 useSolve,
33+ useWebGpuSupport,
3234 type ParameterSection,
3335 } from "@plate-lab/components";
3436 import { ResultsView } from "./views/ResultsView";
@@ -118,11 +120,13 @@ const isMode = (v: unknown): v is Mode =>
118120 typeof v === "string" && (v === "compare" || METHODS.some((m) => m.id === v));
119121 const isPair = (v: unknown): v is [string, string] =>
120122 Array.isArray(v) && v.length === 2 && v.every((x) => METHODS.some((m) => m.id === x));
121-/** Older stored specs lack sections added later; fill them from the defaults. */
122-const migrateSpec = (v: unknown): ProblemSpec => ({
123- ...defaultSpec(),
124- ...(v as Partial<ProblemSpec>),
125-});
123+/** Older stored specs lack sections and fields added later; fill them from the defaults. */
124+const migrateSpec = (v: unknown): ProblemSpec => {
125+ const d = defaultSpec();
126+ const s = (v ?? {}) as Partial<ProblemSpec>;
127+ // plate.shape arrived after the first release; a stored plate without it is a rectangle
128+ return { ...d, ...s, plate: { ...d.plate, ...(s.plate ?? {}) } };
129+};
126130 const isRightTab = (v: unknown): v is RightTab => RIGHT_TABS.some((t) => t.id === v);
127131 const isBool = (v: unknown): v is boolean => typeof v === "boolean";
128132 const isEdits = (v: unknown): v is Record<string, string> =>
@@ -154,6 +158,10 @@ export function App() {
154158 () => false,
155159 isBool,
156160 );
161+ // The explicit scheme's loop can run as its WebGPU port; off by default, the numbl loop being the reference.
162+ const [fdGpu, setFdGpu] = usePersistentState<boolean>("fdGpu", () => false, isBool);
163+ const gpuSupport = useWebGpuSupport();
164+ const gpuOn = fdGpu && gpuSupport.status === "available";
157165 const runner = useRunner();
158166 const [backend, setBackend] = useState<BackendInfo | null>(null);
159167 useEffect(() => {
@@ -185,12 +193,19 @@ export function App() {
185193 [method, methodA, methodB],
186194 );
187195
188- // Source bundles: what the code view shows is what runs
196+ const mainEdited = ["main.m", "main_a.m", "main_b.m"].some((n) => n in edits);
197+
198+ // Source bundles: what the code view shows is what runs. With the GPU option on, a method
199+ // with a ported loop gets its script written as the three stages the host then runs, and
200+ // the stage scripts themselves; an edited main script runs as written, in numbl.
189201 const bundles = useMemo(() => {
202+ const gpuFor = (m: MethodDefinition) => gpuOn && !!m.gpu;
203+ const stagesFor = (s: ProblemSpec, m: MethodDefinition) =>
204+ gpuFor(m) && !mainEdited ? (buildGpuStages(s, m) ?? undefined) : undefined;
190205 if (mode === "compare") {
191206 const [specA, specB] = matchedSpecs(spec, methodA, methodB);
192- const a = renameEntry(buildRunBundle(specA, methodA), "main_a.m");
193- const b = renameEntry(buildRunBundle(specB, methodB), "main_b.m");
207+ const a = renameEntry(buildRunBundle(specA, methodA, { gpu: gpuFor(methodA) }), "main_a.m");
208+ const b = renameEntry(buildRunBundle(specB, methodB, { gpu: gpuFor(methodB) }), "main_b.m");
194209 const c = buildCompareBundle();
195210 const display: SourceBundle = {
196211 files: [a.files[0], b.files[0], c.files[0], ...LIBRARY_FILES],
@@ -202,13 +217,13 @@ export function App() {
202217 runA: withEdits(a, edits),
203218 runB: withEdits(b, edits),
204219 runC: withEdits(c, edits),
220+ gpuA: stagesFor(specA, methodA),
221+ gpuB: stagesFor(specB, methodB),
205222 };
206223 }
207- const single = buildRunBundle(spec, method!);
208- return { display: single, run: withEdits(single, edits) };
209- }, [mode, spec, edits, method, methodA, methodB]);
210-
211- const mainEdited = ["main.m", "main_a.m", "main_b.m"].some((n) => n in edits);
224+ const single = buildRunBundle(spec, method!, { gpu: gpuFor(method!) });
225+ return { display: single, run: withEdits(single, edits), gpu: stagesFor(spec, method!) };
226+ }, [mode, spec, edits, method, methodA, methodB, gpuOn, mainEdited]);
212227
213228 // Two result slots: the single method uses A; compare uses A and B
214229 const solveA = useSolve(parsePlateResult);
@@ -236,17 +251,22 @@ export function App() {
236251
237252 const observe = (
238253 m: MethodDefinition,
239- r: { ok: boolean; wallMs: number; value?: unknown } | undefined,
254+ r: { ok: boolean; wallMs: number; value?: unknown; compute?: { loop: string } } | undefined,
240255 ) => {
241256 if (!r?.ok || !m.planFunction) return;
242257 const res = parsePlateResult(r.value as never);
243- cost.current.observe(m, Number(res.info.extra.cell_steps), r.wallMs);
258+ cost.current.observe(
259+ m,
260+ Number(res.info.extra.cell_steps),
261+ r.wallMs,
262+ r.compute?.loop === "webgpu",
263+ );
244264 };
245265
246266 /** Identifies the inputs a result came from, so a result from other settings can say so. */
247267 const stamp = useMemo(
248- () => JSON.stringify({ spec, edits, mode, pair }),
249- [spec, edits, mode, pair],
268+ () => JSON.stringify({ spec, edits, mode, pair, gpuOn }),
269+ [spec, edits, mode, pair, gpuOn],
250270 );
251271 const [ranWith, setRanWith] = useState<{ mode: Mode; stamp: string } | null>(null);
252272
@@ -259,13 +279,13 @@ export function App() {
259279 solveCmp.reset();
260280 setRanWith({ mode, stamp });
261281 if (method) {
262- const r = await solveA.run({ bundle: bundles.run! });
282+ const r = await solveA.run({ bundle: bundles.run!, gpu: bundles.gpu });
263283 observe(method, r);
264284 } else {
265- const a = await solveA.run({ bundle: bundles.runA! });
285+ const a = await solveA.run({ bundle: bundles.runA!, gpu: bundles.gpuA });
266286 if (!a?.ok) return;
267287 observe(methodA, a);
268- const b = await solveB.run({ bundle: bundles.runB! });
288+ const b = await solveB.run({ bundle: bundles.runB!, gpu: bundles.gpuB });
269289 if (!b?.ok) return;
270290 observe(methodB, b);
271291 await solveCmp.run({ bundle: bundles.runC!, inputs: { A: a.value, B: b.value } });
@@ -300,6 +320,7 @@ export function App() {
300320 setEdits({});
301321 setRightTab("equations");
302322 setAdvancedOpen(false);
323+ setFdGpu(false);
303324 solveA.reset();
304325 solveB.reset();
305326 solveCmp.reset();
@@ -339,7 +360,7 @@ export function App() {
339360 (currentResult !== null || solveCmp.last !== null);
340361 const estimateMs =
341362 planMethod && methodPlan
342- ? cost.current.estimateMs(planMethod, methodPlan.cell_steps)
363+ ? cost.current.estimateMs(planMethod, methodPlan.cell_steps, gpuOn && !!planMethod.gpu)
343364 : undefined;
344365 const slow = estimateMs !== undefined && estimateMs > 5000;
345366 const advancedIssue = issues.some(
@@ -527,12 +548,13 @@ export function App() {
527548 fdDtMax={planMethod?.id === "fd" ? methodPlan?.dt_max : undefined}
528549 sections={methodSections}
529550 open={methodSections}
551+ fdGpu={{ enabled: fdGpu, onChange: setFdGpu, support: gpuSupport }}
530552 />
531553 {methodPlan && planMethod && (
532554 <div className="pl-small pl-muted plan-line">
533555 {planMethod.shortName}: {methodPlan.nsteps.toLocaleString()} steps of{" "}
534556 {(methodPlan.dt * 1e6).toPrecision(3)} µs, {methodPlan.cell_steps.toLocaleString()}{" "}
535- cell-steps.
557+ cell-steps{gpuOn && planMethod.gpu ? ", time loop on the GPU" : ""}.
536558 </div>
537559 )}
538560 {/* The modal method has no grid or step of its own, so its simulated window lives here. */}
@@ -778,6 +800,7 @@ function CompareResults({
778800 frameIndex={idx}
779801 Lx={spec.plate.Lx}
780802 Ly={spec.plate.Ly}
803+ shape={spec.plate.shape}
781804 title={`w, ${la}`}
782805 unit="µm"
783806 scale={1e6}
@@ -789,6 +812,7 @@ function CompareResults({
789812 frameIndex={idx}
790813 Lx={spec.plate.Lx}
791814 Ly={spec.plate.Ly}
815+ shape={spec.plate.shape}
792816 title={`w, ${lb}`}
793817 unit="µm"
794818 scale={1e6}
@@ -802,6 +826,7 @@ function CompareResults({
802826 frameIndex={idx}
803827 Lx={spec.plate.Lx}
804828 Ly={spec.plate.Ly}
829+ shape={spec.plate.shape}
805830 title={`difference, ${la} − ${lb}`}
806831 unit="µm"
807832 scale={1e6}
@@ -816,6 +841,7 @@ function CompareResults({
816841 frameIndex={idx}
817842 Lx={spec.plate.Lx}
818843 Ly={spec.plate.Ly}
844+ shape={spec.plate.shape}
819845 title={`von Mises, ${la}`}
820846 unit="MPa"
821847 scale={1e-6}
@@ -826,6 +852,7 @@ function CompareResults({
826852 frameIndex={idx}
827853 Lx={spec.plate.Lx}
828854 Ly={spec.plate.Ly}
855+ shape={spec.plate.shape}
829856 title={`von Mises, ${lb}`}
830857 unit="MPa"
831858 scale={1e-6}
apps/lab/src/views/ResultsView.tsxmodified+24−6View file
@@ -75,13 +75,18 @@ export function ResultsView({ result, spec, method, fieldWidth = 460 }: ResultsV
7575 : undefined;
7676 const tFrame = result.frames.t[Math.min(player.index, nf - 1)] ?? 0;
7777
78- const markers: FieldMarker[] = useMemo(
79- () => [
78+ const supportX = result.info.extra.support_x;
79+ const supportY = result.info.extra.support_y;
80+ const markers: FieldMarker[] = useMemo(() => {
81+ const m: FieldMarker[] = [
8082 { x: spec.excitation.x, y: spec.excitation.y, label: "strike", color: "#d1495b" },
8183 { x: result.probe.x, y: result.probe.y, label: "probe", color: "#1f5fa8" },
82- ],
83- [spec.excitation.x, spec.excitation.y, result.probe.x, result.probe.y],
84- );
84+ ];
85+ // a centre-pinned plate reports the node it is held at
86+ if (typeof supportX === "number" && typeof supportY === "number")
87+ m.push({ x: supportX, y: supportY, label: "support", color: "#2a9d8f" });
88+ return m;
89+ }, [spec.excitation.x, spec.excitation.y, result.probe.x, result.probe.y, supportX, supportY]);
8590
8691 const onPlotClick = (t: number) => {
8792 let best = 0;
@@ -102,7 +107,12 @@ export function ResultsView({ result, spec, method, fieldWidth = 460 }: ResultsV
102107 return (
103108 <div className="results">
104109 <div className="results-toolbar sound-bar">
105- <AudioPlayer p={result.audio.p} fs={result.audio.fs} label="sound" />
110+ <AudioPlayer
111+ p={result.audio.p}
112+ fs={result.audio.fs}
113+ label="sound"
114+ filename={`plate-${result.method}`}
115+ />
106116 </div>
107117 <div className="results-toolbar">
108118 <label className="pl-small pl-muted">show</label>
@@ -144,6 +154,7 @@ export function ResultsView({ result, spec, method, fieldWidth = 460 }: ResultsV
144154 frameIndex={player.index}
145155 Lx={spec.plate.Lx}
146156 Ly={spec.plate.Ly}
157+ shape={spec.plate.shape}
147158 title={q.label}
148159 unit={q.unit}
149160 scale={q.scale}
@@ -159,6 +170,13 @@ export function ResultsView({ result, spec, method, fieldWidth = 460 }: ResultsV
159170 : "Surface stress at z = h/2 from the curvature of w; the opposite face carries the opposite sign."
160171 }
161172 />
173+ {spec.plate.shape === "ellipse" && result.method !== "modal" && (
174+ <div className="pl-small pl-muted">
175+ The solver's grid spans the box around the ellipse; the edge is a staircase of grid nodes,
176+ drawn here clipped to the true outline. Stress on the outermost ring is copied from the
177+ nearest node with a full stencil.
178+ </div>
179+ )}
162180 {frames && <FramePlayer {...player} times={result.frames.t} />}
163181
164182 <LinePlot
package-lock.jsonmodified+8−0View file
@@ -16,6 +16,7 @@
1616 "@types/react": "^19.2.5",
1717 "@types/react-dom": "^19.2.3",
1818 "@vitejs/plugin-react": "^5.1.1",
19+ "@webgpu/types": "^0.1.72",
1920 "prettier": "^3.8.1",
2021 "typescript": "~5.9.3",
2122 "vite": "^7.2.4",
@@ -1799,6 +1800,13 @@
17991800 "url": "https://opencollective.com/vitest"
18001801 }
18011802 },
1803+ "node_modules/@webgpu/types": {
1804+ "version": "0.1.72",
1805+ "resolved": "https://registry.npmjs.org/@webgpu/types/-/types-0.1.72.tgz",
1806+ "integrity": "sha512-0cF7RFM2edNoiIS1ODJp0/Gzv4/xSXhwoR0YCza+OWpJWtn4wmo9DvK91aLlH9+uUnwIriP7ZiC3WitmyhuzBw==",
1807+ "dev": true,
1808+ "license": "BSD-3-Clause"
1809+ },
18021810 "node_modules/acorn": {
18031811 "version": "8.18.0",
18041812 "resolved": "https://registry.npmjs.org/acorn/-/acorn-8.18.0.tgz",
package.jsonmodified+1−0View file
@@ -21,6 +21,7 @@
2121 "@types/react": "^19.2.5",
2222 "@types/react-dom": "^19.2.3",
2323 "@vitejs/plugin-react": "^5.1.1",
24+ "@webgpu/types": "^0.1.72",
2425 "prettier": "^3.8.1",
2526 "typescript": "~5.9.3",
2627 "vite": "^7.2.4",
packages/components/src/audio/AudioPlayer.tsxmodified+34−6View file
@@ -1,19 +1,22 @@
11 import { useEffect, useMemo, useRef, useState } from "react";
2-import { prepareAudio, toAudioBuffer } from "@plate-lab/core";
2+import { encodeWav, prepareAudio, toAudioBuffer } from "@plate-lab/core";
33
44 export interface AudioPlayerProps {
55 /** Pressure signal [Pa] */
66 p: Float64Array;
77 fs: number;
88 label?: string;
9+ /** Base name for the exported .wav file (no extension) */
10+ filename?: string;
911 }
1012
1113 /**
12- * Plays the radiated-pressure estimate once, exactly as long as it is. Playback
13- * level is normalized; the displayed peak is not. A clip of a few milliseconds
14+ * Plays the radiated-pressure estimate once, exactly as long as it is, and
15+ * offers the same clip as a 16-bit mono WAV download. Playback and export level
16+ * are normalized; the displayed peak is not. A clip of a few milliseconds
1417 * sounds like a click, and is labelled as such rather than repeated.
1518 */
16-export function AudioPlayer({ p, fs, label }: AudioPlayerProps) {
19+export function AudioPlayer({ p, fs, label, filename = "plate-sound" }: AudioPlayerProps) {
1720 const prep = useMemo(() => prepareAudio(p, fs), [p, fs]);
1821 const ctxRef = useRef<AudioContext | null>(null);
1922 const srcRef = useRef<AudioBufferSourceNode | null>(null);
@@ -38,15 +41,40 @@ export function AudioPlayer({ p, fs, label }: AudioPlayerProps) {
3841 srcRef.current = src;
3942 setPlaying(true);
4043 };
44+ const download = () => {
45+ // The exported file carries the same normalized signal that is played, at the
46+ // simulation's own sample rate; the peak in Pa is reported alongside it.
47+ const blob = new Blob([encodeWav(prep.samples, Math.round(prep.fs)) as BlobPart], {
48+ type: "audio/wav",
49+ });
50+ const url = URL.createObjectURL(blob);
51+ const a = document.createElement("a");
52+ a.href = url;
53+ a.download = `${filename}.wav`;
54+ a.click();
55+ setTimeout(() => URL.revokeObjectURL(url), 0);
56+ };
57+
4158 const short = prep.durationS < 0.1;
4259 return (
4360 <div className="pl-row pl-small">
4461 <button className="pl-btn" onClick={playing ? stop : play} disabled={p.length < 2}>
4562 {playing ? "Stop" : `Play${label ? ` ${label}` : ""}`}
4663 </button>
64+ <button
65+ className="pl-btn"
66+ onClick={download}
67+ disabled={p.length < 2}
68+ title="Download as a 16-bit mono WAV file"
69+ >
70+ Save .wav
71+ </button>
4772 <span className="pl-muted">
48- {(prep.durationS * 1000).toFixed(1)} ms at {fs} Hz, peak {prep.peakPa.toExponential(2)} Pa; normalized for playback, played once
49- {short ? ". This excerpt is too short to hear as a note; lengthen the simulated window for a longer clip" : ""}
73+ {(prep.durationS * 1000).toFixed(1)} ms at {fs} Hz, peak {prep.peakPa.toExponential(2)} Pa;
74+ normalized for playback, played once
75+ {short
76+ ? ". This excerpt is too short to hear as a note; lengthen the simulated window for a longer clip"
77+ : ""}
5078 </span>
5179 </div>
5280 );
packages/components/src/field/FieldView.tsxmodified+58−8View file
@@ -7,6 +7,7 @@ import {
77 symmetricRange,
88 type Field2D,
99 type FieldFrames,
10+ type PlateShape,
1011 type Range,
1112 } from "@plate-lab/core";
1213 import { Colorbar, fmt } from "./Colorbar";
@@ -23,9 +24,11 @@ export interface FieldViewProps {
2324 field?: Field2D;
2425 frames?: FieldFrames;
2526 frameIndex?: number;
26- /** Physical extents [m] */
27+ /** Physical extents of the bounding box [m] */
2728 Lx: number;
2829 Ly: number;
30+ /** Outline of the plate within the box; default rectangle */
31+ shape?: PlateShape;
2932 title?: string;
3033 unit?: string;
3134 /** Display multiplier applied to values in labels (e.g. 1e-6 for MPa) */
@@ -43,12 +46,30 @@ export interface FieldViewProps {
4346
4447 /**
4548 * Paints a field on the plate outline. Values pass straight from the solver
46- * through a colormap; nothing here changes them.
49+ * through a colormap; nothing here changes them. Nodes the solver marked NaN
50+ * (outside an elliptical plate) come out transparent, and the image is clipped
51+ * to the plate's true outline.
4752 */
4853 export function FieldView(props: FieldViewProps) {
49- const { frames, frameIndex = 0, Lx, Ly, title, unit, scale = 1, colormap, symmetric, markers = [], width = 420, caption } = props;
50- const field: Field2D | undefined = props.field ?? (frames ? frameOf(frames, Math.min(frameIndex, frames.nf - 1)) : undefined);
54+ const {
55+ frames,
56+ frameIndex = 0,
57+ Lx,
58+ Ly,
59+ shape = "rectangle",
60+ title,
61+ unit,
62+ scale = 1,
63+ colormap,
64+ symmetric,
65+ markers = [],
66+ width = 420,
67+ caption,
68+ } = props;
69+ const field: Field2D | undefined =
70+ props.field ?? (frames ? frameOf(frames, Math.min(frameIndex, frames.nf - 1)) : undefined);
5171 const cmId = colormap ?? (symmetric ? "bwr" : "viridis");
72+ const ellipse = shape === "ellipse";
5273
5374 const globalRange = useMemo(() => {
5475 const r = frames ? rangeOf(frames.data) : field ? rangeOf(field.data) : { min: 0, max: 1 };
@@ -86,18 +107,37 @@ export function FieldView(props: FieldViewProps) {
86107 pixRef.current = pix;
87108 const octx = off.getContext("2d");
88109 if (!octx) return;
89- octx.putImageData(new ImageData(pix as Uint8ClampedArray<ArrayBuffer>, field.nx, field.ny), 0, 0);
110+ octx.putImageData(
111+ new ImageData(pix as Uint8ClampedArray<ArrayBuffer>, field.nx, field.ny),
112+ 0,
113+ 0,
114+ );
90115 ctx.save();
91116 ctx.scale(dpr, dpr);
92117 ctx.clearRect(0, 0, width, height);
93118 ctx.imageSmoothingEnabled = true;
119+ const outline = () => {
120+ ctx.beginPath();
121+ ctx.ellipse(width / 2, height / 2, width / 2 - 0.75, height / 2 - 0.75, 0, 0, 2 * Math.PI);
122+ };
123+ if (ellipse) {
124+ ctx.save();
125+ outline();
126+ ctx.clip();
127+ }
94128 // Grid nodes sit at cell centres of the image; map node (0..nx-1) to the full extent
95129 const sx = width / (field.nx - 1);
96130 const sy = height / (field.ny - 1);
97131 ctx.drawImage(off, -sx / 2, -sy / 2, width + sx, height + sy);
132+ if (ellipse) ctx.restore();
98133 ctx.strokeStyle = "#333";
99134 ctx.lineWidth = 1;
100- ctx.strokeRect(0.5, 0.5, width - 1, height - 1);
135+ if (ellipse) {
136+ outline();
137+ ctx.stroke();
138+ } else {
139+ ctx.strokeRect(0.5, 0.5, width - 1, height - 1);
140+ }
101141 for (const m of markers) {
102142 const px = (m.x / Lx) * width;
103143 const py = height - (m.y / Ly) * height;
@@ -116,7 +156,7 @@ export function FieldView(props: FieldViewProps) {
116156 ctx.fillText(m.label, px + 8, py - 6);
117157 }
118158 ctx.restore();
119- }, [field, range, cmId, width, height, markers, Lx, Ly]);
159+ }, [field, range, cmId, width, height, markers, Lx, Ly, ellipse]);
120160
121161 const peak = field ? rangeOf(field.data) : null;
122162
@@ -131,7 +171,17 @@ export function FieldView(props: FieldViewProps) {
131171 )}
132172 </div>
133173 <div className="pl-fieldview-body">
134- <canvas ref={canvasRef} style={{ width, height, aspectRatio: `${Lx} / ${Ly}` }} aria-label={title} />
174+ <canvas
175+ ref={canvasRef}
176+ style={{
177+ width,
178+ height,
179+ aspectRatio: `${Lx} / ${Ly}`,
180+ border: ellipse ? "none" : undefined,
181+ background: ellipse ? "transparent" : undefined,
182+ }}
183+ aria-label={title}
184+ />
135185 <Colorbar range={range} colormap={cmId} unit={unit} scale={scale} height={height} />
136186 </div>
137187 {caption && <div className="pl-small pl-muted">{caption}</div>}
packages/components/src/hooks/useWebGpu.tsadded+39−0View file
@@ -0,0 +1,39 @@
1+import { useEffect, useState } from "react";
2+import { describeAdapter } from "@plate-lab/core";
3+
4+export type WebGpuSupport =
5+ | { status: "checking" }
6+ | { status: "available"; adapter: string }
7+ | { status: "unavailable"; reason: string };
8+
9+/**
10+ * Whether this browser exposes WebGPU, and on what. Asked once on the main
11+ * thread so the controls can say what a GPU run would use; the worker makes
12+ * its own device request when a run asks for it.
13+ */
14+export function useWebGpuSupport(): WebGpuSupport {
15+ const [support, setSupport] = useState<WebGpuSupport>({ status: "checking" });
16+ useEffect(() => {
17+ let alive = true;
18+ const set = (s: WebGpuSupport) => {
19+ if (alive) setSupport(s);
20+ };
21+ (async () => {
22+ try {
23+ const gpu = (navigator as { gpu?: GPU }).gpu;
24+ if (!gpu)
25+ return set({ status: "unavailable", reason: "this browser does not expose WebGPU" });
26+ const adapter = await gpu.requestAdapter();
27+ if (!adapter)
28+ return set({ status: "unavailable", reason: "no WebGPU adapter was offered" });
29+ set({ status: "available", adapter: describeAdapter(adapter) });
30+ } catch (e) {
31+ set({ status: "unavailable", reason: e instanceof Error ? e.message : String(e) });
32+ }
33+ })();
34+ return () => {
35+ alive = false;
36+ };
37+ }, []);
38+ return support;
39+}
packages/components/src/index.tsmodified+1−0View file
@@ -1,5 +1,6 @@
11 export { RunnerProvider, useRunner } from "./hooks/RunnerContext";
22 export { useSolve, type Solve, type SolveState } from "./hooks/useSolve";
3+export { useWebGpuSupport, type WebGpuSupport } from "./hooks/useWebGpu";
34 export { CodeView, type CodeViewProps } from "./code/CodeView";
45 export { NumberField, type NumberFieldProps } from "./params/NumberField";
56 export { CflSlider, type CflSliderProps } from "./params/CflSlider";
packages/components/src/params/ParameterPanel.tsxmodified+92−6View file
@@ -5,11 +5,13 @@ import {
55 type BoundaryCondition,
66 type Excitation,
77 type MethodDefinition,
8+ type PlateShape,
89 type ProblemSpec,
910 type ValidationIssue,
1011 } from "@plate-lab/core";
1112 import { NumberField } from "./NumberField";
1213 import { CflSlider } from "./CflSlider";
14+import type { WebGpuSupport } from "../hooks/useWebGpu";
1315
1416 export type ParameterSection =
1517 | "plate"
@@ -50,11 +52,22 @@ export interface ParameterPanelProps {
5052 fdDtMax?: number;
5153 /** Sections open initially; default all */
5254 open?: ParameterSection[];
55+ /** The GPU option of the explicit scheme: whether it is on, and what this browser offers */
56+ fdGpu?: { enabled: boolean; onChange: (v: boolean) => void; support: WebGpuSupport };
5357 }
5458
5559 /** Controls that edit the spec. Every control maps to one field of the `params` struct. */
5660 export function ParameterPanel(props: ParameterPanelProps) {
57- const { spec, onChange, sections = ALL_SECTIONS, method, issues = [], fdDtMax, open } = props;
61+ const {
62+ spec,
63+ onChange,
64+ sections = ALL_SECTIONS,
65+ method,
66+ issues = [],
67+ fdDtMax,
68+ open,
69+ fdGpu,
70+ } = props;
5871 const set = (patch: Partial<ProblemSpec>) => onChange({ ...spec, ...patch });
5972 const issue = (path: string) =>
6073 issues.find((i) => i.path === path || i.path.startsWith(path + "."));
@@ -68,8 +81,12 @@ export function ParameterPanel(props: ParameterPanelProps) {
6881 )?.id;
6982 const platePresetId = PLATE_PRESETS.find(
7083 (p) =>
71- p.plate.Lx === spec.plate.Lx && p.plate.Ly === spec.plate.Ly && p.plate.h === spec.plate.h,
84+ p.plate.shape === spec.plate.shape &&
85+ p.plate.Lx === spec.plate.Lx &&
86+ p.plate.Ly === spec.plate.Ly &&
87+ p.plate.h === spec.plate.h,
7288 )?.id;
89+ const ellipse = spec.plate.shape === "ellipse";
7390
7491 return (
7592 <div className="pl-params">
@@ -98,8 +115,27 @@ export function ParameterPanel(props: ParameterPanelProps) {
98115 {PLATE_PRESETS.find((p) => p.id === platePresetId)?.note}
99116 </div>
100117 )}
118+ <div className="pl-field pl-field-wide">
119+ <label>Shape</label>
120+ <select
121+ className="pl-select"
122+ value={spec.plate.shape}
123+ onChange={(ev) =>
124+ set({ plate: { ...spec.plate, shape: ev.target.value as PlateShape } })
125+ }
126+ >
127+ <option value="rectangle">rectangle</option>
128+ <option value="ellipse">ellipse (a disk when Lx = Ly)</option>
129+ </select>
130+ </div>
131+ {method && !method.supports.shape.includes(spec.plate.shape) && (
132+ <div className="pl-issue">
133+ {method.name} has no {spec.plate.shape} solution; the finite-difference methods embed
134+ it in their grid.
135+ </div>
136+ )}
101137 <NumberField
102- label="Length Lx"
138+ label={ellipse ? "Diameter along x, Lx" : "Length Lx"}
103139 value={spec.plate.Lx}
104140 scale={1000}
105141 unit="mm"
@@ -107,7 +143,7 @@ export function ParameterPanel(props: ParameterPanelProps) {
107143 onChange={(v) => set({ plate: { ...spec.plate, Lx: v } })}
108144 />
109145 <NumberField
110- label="Width Ly"
146+ label={ellipse ? "Diameter along y, Ly" : "Width Ly"}
111147 value={spec.plate.Ly}
112148 scale={1000}
113149 unit="mm"
@@ -123,6 +159,13 @@ export function ParameterPanel(props: ParameterPanelProps) {
123159 invalid={invalid("plate.h")}
124160 onChange={(v) => set({ plate: { ...spec.plate, h: v } })}
125161 />
162+ {ellipse && (
163+ <div className="pl-small pl-muted">
164+ The ellipse fills the Lx by Ly box, and positions are still measured from the box's
165+ corner. The finite-difference grid spans the box; the edge becomes a staircase of
166+ nodes, first order in the grid spacing.
167+ </div>
168+ )}
126169 </Group>
127170 )}
128171 {show("material") && (
@@ -193,6 +236,9 @@ export function ParameterPanel(props: ParameterPanelProps) {
193236 <option value="simply-supported">simply supported (w = 0, no moment)</option>
194237 <option value="clamped">clamped (w = 0, no slope)</option>
195238 <option value="free">free (no moment, no shear): the free-free plate</option>
239+ <option value="pinned-center">
240+ free edges, pinned at the centre: a cymbal on its stand
241+ </option>
196242 </select>
197243 </div>
198244 {method && !method.supports.bc.includes(spec.bc) && (
@@ -202,8 +248,8 @@ export function ParameterPanel(props: ParameterPanelProps) {
202248 )}
203249 <div className="pl-small pl-muted">
204250 All three are idealizations: a real clamp has finite stiffness, a real support some
205- rotational restraint, and a free plate rests on something that takes its rigid-body
206- motion, which the solvers remove.
251+ rotational restraint, a free plate rests on something that takes its rigid-body motion,
252+ which the solvers remove, and a stand holds a cymbal over a washer, not at a point.
207253 </div>
208254 </Group>
209255 )}
@@ -473,6 +519,7 @@ export function ParameterPanel(props: ParameterPanelProps) {
473519 dtMax={fdDtMax}
474520 onChange={(v) => set({ fd: { ...spec.fd, cfl: v } })}
475521 />
522+ {fdGpu && <GpuToggle {...fdGpu} />}
476523 </Group>
477524 )}
478525 {show("implicit") && (
@@ -534,6 +581,45 @@ export function ParameterPanel(props: ParameterPanelProps) {
534581 );
535582 }
536583
584+/**
585+ * The explicit scheme's time loop has a WebGPU port. The option is off by
586+ * default: the double-precision loop in numbl is the reference, and the port is
587+ * the same arithmetic in single precision, many times faster.
588+ */
589+function GpuToggle({
590+ enabled,
591+ onChange,
592+ support,
593+}: {
594+ enabled: boolean;
595+ onChange: (v: boolean) => void;
596+ support: WebGpuSupport;
597+}) {
598+ const available = support.status === "available";
599+ return (
600+ <div className="pl-gpu">
601+ <label className={`pl-check${available ? "" : " pl-disabled"}`}>
602+ <input
603+ type="checkbox"
604+ checked={enabled && available}
605+ disabled={!available}
606+ onChange={(ev) => onChange(ev.target.checked)}
607+ />
608+ Run the time loop on the GPU (WebGPU)
609+ </label>
610+ <div className="pl-small pl-muted">
611+ {support.status === "checking" && "Looking for WebGPU…"}
612+ {support.status === "unavailable" &&
613+ `Not available here (${support.reason}); the loop runs in numbl.`}
614+ {support.status === "available" &&
615+ (enabled
616+ ? `plate_fd_loop.m runs as its WebGPU port on ${support.adapter}, in single precision; setup, stresses, sound and spectra stay in numbl. Editing plate_fd_loop.m sends the loop back to numbl.`
617+ : `Available on ${support.adapter}. The loop then runs in single precision, typically 10 to 50 times faster than in numbl.`)}
618+ </div>
619+ </div>
620+ );
621+}
622+
537623 function customName(name: string): string {
538624 return name.endsWith(" (edited)") ? name : `${name} (edited)`;
539625 }
packages/components/src/status/RunStatus.tsxmodified+60−16View file
@@ -12,9 +12,19 @@ export interface RunStatusProps {
1212 }
1313
1414 /** Makes the execution visible: what ran, where, how big it was and how long it took. */
15-export function RunStatus({ state, result, method, backend, estimateMs, onCancel }: RunStatusProps) {
15+export function RunStatus({
16+ state,
17+ result,
18+ method,
19+ backend,
20+ estimateMs,
21+ onCancel,
22+}: RunStatusProps) {
1623 const be = state.status === "done" || state.status === "error" ? state.run.backend : backend;
1724 const wallMs = state.status === "done" || state.status === "error" ? state.run.wallMs : null;
25+ const compute =
26+ state.status === "done" || state.status === "error" ? state.run.compute : undefined;
27+ const onGpu = compute?.loop === "webgpu";
1828 return (
1929 <div className="pl-status">
2030 <div className="pl-row">
@@ -30,8 +40,16 @@ export function RunStatus({ state, result, method, backend, estimateMs, onCancel
3040 </>
3141 )}
3242 {state.status === "done" && <span className="pl-badge pl-badge-ok">finished</span>}
33- {state.status === "error" && <span className="pl-badge pl-badge-error">{state.run.cancelled ? "cancelled" : "error"}</span>}
34- {result && !result.info.stable && <span className="pl-badge pl-badge-error">diverged at step {result.info.diverged_step}</span>}
43+ {state.status === "error" && (
44+ <span className="pl-badge pl-badge-error">
45+ {state.run.cancelled ? "cancelled" : "error"}
46+ </span>
47+ )}
48+ {result && !result.info.stable && (
49+ <span className="pl-badge pl-badge-error">
50+ diverged at step {result.info.diverged_step}
51+ </span>
52+ )}
3553 {method && <span className="pl-muted">{method.name}</span>}
3654 </div>
3755 {state.status === "running" && (
@@ -45,17 +63,28 @@ export function RunStatus({ state, result, method, backend, estimateMs, onCancel
4563 <tr>
4664 <td>Backend</td>
4765 <td>
48- {be ? `${be.engine} ${be.version}, ${be.optimization === "1" ? "JS JIT (opt 1)" : "interpreter (opt 0)"}, ${be.host === "web-worker" ? "Web Worker in this browser" : "Node.js"}` : "numbl (not yet started)"}
66+ {be
67+ ? `${be.engine} ${be.version}, ${be.optimization === "1" ? "JS JIT (opt 1)" : "interpreter (opt 0)"}, ${be.host === "web-worker" ? "Web Worker in this browser" : "Node.js"}`
68+ : "numbl (not yet started)"}
4969 </td>
5070 </tr>
5171 <tr>
5272 <td>Wall time</td>
5373 <td>
5474 {wallMs !== null ? (
55- <>
56- {fmtMs(wallMs)} in numbl
57- {result ? ` (${fmtMs(result.info.wall_s * 1000)} inside the solver by tic/toc)` : ""}
58- </>
75+ onGpu ? (
76+ <>
77+ {fmtMs(wallMs)} in total: {fmtMs(compute.gpuMs ?? 0)} for the time loop on the
78+ GPU, {fmtMs(compute.numblMs ?? 0)} in numbl for setup and output
79+ </>
80+ ) : (
81+ <>
82+ {fmtMs(wallMs)} in numbl
83+ {result
84+ ? ` (${fmtMs(result.info.wall_s * 1000)} inside the solver by tic/toc)`
85+ : ""}
86+ </>
87+ )
5988 ) : state.status === "running" ? (
6089 "…"
6190 ) : estimateMs !== undefined ? (
@@ -65,6 +94,16 @@ export function RunStatus({ state, result, method, backend, estimateMs, onCancel
6594 )}
6695 </td>
6796 </tr>
97+ {compute && (
98+ <tr>
99+ <td>Time loop</td>
100+ <td>
101+ {onGpu
102+ ? `WebGPU port of plate_fd_loop.m on ${compute.adapter ?? "the GPU"}, single precision`
103+ : `numbl${compute.note ? `: ${compute.note}` : ""}`}
104+ </td>
105+ </tr>
106+ )}
68107 {result && (
69108 <tr>
70109 <td>Problem size</td>
@@ -75,7 +114,8 @@ export function RunStatus({ state, result, method, backend, estimateMs, onCancel
75114 <tr>
76115 <td>Time step</td>
77116 <td>
78- {fmtTime(result.info.dt)}, {result.info.nsteps.toLocaleString()} steps, {result.t.length.toLocaleString()} samples at the probe
117+ {fmtTime(result.info.dt)}, {result.info.nsteps.toLocaleString()} steps,{" "}
118+ {result.t.length.toLocaleString()} samples at the probe
79119 </td>
80120 </tr>
81121 )}
@@ -90,19 +130,23 @@ export function RunStatus({ state, result, method, backend, estimateMs, onCancel
90130 ))}
91131 </div>
92132 )}
93- {(state.status === "running" || state.status === "done" || state.status === "error") && state.output.length > 0 && (
94- <details>
95- <summary className="pl-small pl-muted">Script output ({state.output.length} lines)</summary>
96- <div className="pl-log">{state.output.join("")}</div>
97- </details>
98- )}
133+ {(state.status === "running" || state.status === "done" || state.status === "error") &&
134+ state.output.length > 0 && (
135+ <details>
136+ <summary className="pl-small pl-muted">
137+ Script output ({state.output.length} lines)
138+ </summary>
139+ <div className="pl-log">{state.output.join("")}</div>
140+ </details>
141+ )}
99142 </div>
100143 );
101144 }
102145
103146 function problemSize(r: PlateResult): string {
104147 const x = r.info.extra;
105- if (r.method === "modal") return `${x.nmodes} modes, ${r.info.ncells.toLocaleString()} display points, ${r.frames.t.length} frames`;
148+ if (r.method === "modal")
149+ return `${x.nmodes} modes, ${r.info.ncells.toLocaleString()} display points, ${r.frames.t.length} frames`;
106150 if (r.method === "fd")
107151 return `${x.nx} × ${x.ny} = ${r.info.ncells.toLocaleString()} nodes, ${Number(x.cell_steps).toLocaleString()} cell-steps, ${r.frames.t.length} frames`;
108152 if (r.method === "implicit")
packages/components/src/styles.cssmodified+3−0View file
@@ -48,6 +48,9 @@ button.pl-btn-danger { background: #fff; color: var(--pl-error); border-color: v
4848 .pl-issue { font-size: 11px; color: var(--pl-error); margin: 2px 0 4px; }
4949 .pl-issue-warning { color: var(--pl-warn); }
5050 .pl-disabled { opacity: 0.55; }
51+.pl-gpu { margin-top: 8px; padding-top: 6px; border-top: 1px dashed var(--pl-border); }
52+.pl-check { display: flex; gap: 6px; align-items: center; font-size: 12px; cursor: pointer; }
53+.pl-check input { margin: 0; }
5154
5255 /* Code view */
5356 .pl-code { display: flex; flex-direction: column; border: 1px solid var(--pl-border); border-radius: 6px; overflow: hidden; background: #fff; min-height: 0; }
packages/core/src/audio/wav.tsadded+35−0View file
@@ -0,0 +1,35 @@
1+/**
2+ * Encode a mono signal as a 16-bit PCM RIFF/WAVE file. The input is expected
3+ * to be peak-normalized already (see `prepareAudio`); values outside [-1, 1]
4+ * are clipped rather than wrapped. Non-finite samples are written as silence.
5+ */
6+export function encodeWav(samples: ArrayLike<number>, fs: number): Uint8Array {
7+ const n = samples.length;
8+ const bytesPerSample = 2;
9+ const dataBytes = n * bytesPerSample;
10+ const buf = new ArrayBuffer(44 + dataBytes);
11+ const v = new DataView(buf);
12+ const ascii = (off: number, s: string) => {
13+ for (let i = 0; i < s.length; i++) v.setUint8(off + i, s.charCodeAt(i));
14+ };
15+ ascii(0, "RIFF");
16+ v.setUint32(4, 36 + dataBytes, true);
17+ ascii(8, "WAVE");
18+ ascii(12, "fmt ");
19+ v.setUint32(16, 16, true); // fmt chunk size
20+ v.setUint16(20, 1, true); // PCM
21+ v.setUint16(22, 1, true); // mono
22+ v.setUint32(24, fs, true);
23+ v.setUint32(28, fs * bytesPerSample, true); // byte rate
24+ v.setUint16(32, bytesPerSample, true); // block align
25+ v.setUint16(34, 16, true); // bits per sample
26+ ascii(36, "data");
27+ v.setUint32(40, dataBytes, true);
28+ let off = 44;
29+ for (let i = 0; i < n; i++, off += 2) {
30+ const x = samples[i];
31+ const c = Number.isFinite(x) ? Math.max(-1, Math.min(1, x)) : 0;
32+ v.setInt16(off, Math.round(c < 0 ? c * 32768 : c * 32767), true);
33+ }
34+ return new Uint8Array(buf);
35+}
packages/core/src/cost.tsmodified+17−7View file
@@ -4,25 +4,35 @@ import type { MethodDefinition } from "./methods/types";
44 * Wall-time estimate for grid-based methods. This is a calibration, not
55 * physics: the plan (cell count, steps) comes from the method's plan function
66 * running in numbl; only the seconds-per-cell-step constant lives here, seeded
7- * from the method definition and refined by each finished run.
7+ * from the method definition and refined by each finished run. The GPU path of
8+ * a method is calibrated separately: its cost is nearly independent of the
9+ * cell count at the grids used here, so its per-cell-step rate is a stand-in
10+ * that the measurements correct for the grid in use.
811 */
912 export class CostModel {
1013 private us = new Map<string, number>();
1114 private readonly overheadMs = 150;
1215
13- private rate(method: MethodDefinition): number {
14- return this.us.get(method.id) ?? method.usPerCellStep;
16+ private key(method: MethodDefinition, gpu: boolean): string {
17+ return gpu ? `${method.id}:gpu` : method.id;
1518 }
1619
17- estimateMs(method: MethodDefinition, cellSteps: number): number {
18- return (cellSteps * this.rate(method)) / 1000 + this.overheadMs;
20+ private rate(method: MethodDefinition, gpu: boolean): number {
21+ return (
22+ this.us.get(this.key(method, gpu)) ??
23+ (gpu ? (method.gpu?.usPerCellStep ?? method.usPerCellStep) : method.usPerCellStep)
24+ );
25+ }
26+
27+ estimateMs(method: MethodDefinition, cellSteps: number, gpu = false): number {
28+ return (cellSteps * this.rate(method, gpu)) / 1000 + this.overheadMs;
1929 }
2030
2131 /** Blend a measurement in, weighting recent runs more. */
22- observe(method: MethodDefinition, cellSteps: number, wallMs: number): void {
32+ observe(method: MethodDefinition, cellSteps: number, wallMs: number, gpu = false): void {
2333 if (cellSteps <= 0 || wallMs <= 0) return;
2434 const measured = ((wallMs - this.overheadMs) * 1000) / cellSteps;
2535 if (!Number.isFinite(measured) || measured <= 0) return;
26- this.us.set(method.id, 0.5 * this.rate(method) + 0.5 * measured);
36+ this.us.set(this.key(method, gpu), 0.5 * this.rate(method, gpu) + 0.5 * measured);
2737 }
2838 }
packages/core/src/generate/bundle.tsmodified+71−13View file
@@ -1,5 +1,6 @@
11 import type { ProblemSpec } from "../spec/types";
22 import type { MethodDefinition } from "../methods/types";
3+import type { GpuStages } from "../runner/types";
34 import { generateParamsBlock } from "./params";
45 import { LIBRARY_FILES, type SourceFile } from "./library";
56
@@ -14,30 +15,85 @@ export interface SourceBundle {
1415
1516 export const MAIN_FILE = "main.m";
1617
17-/** The script a method runs: the params block followed by the solver call. */
18-export function generateMain(spec: ProblemSpec, method: MethodDefinition): string {
19- return [
20- generateParamsBlock(spec),
21- "",
22- `% Solve. ${method.solverFunction} and the functions it calls are the library files in the other tabs.`,
23- `result = ${method.solverFunction}(params);`,
24- "",
25- ].join("\n");
18+/**
19+ * The script a method runs: the params block followed by the solver call. With
20+ * `gpu`, for a method whose loop has a WebGPU port, the call is written out as
21+ * its three stages, since that is what the host then runs (see buildGpuStages).
22+ */
23+export function generateMain(
24+ spec: ProblemSpec,
25+ method: MethodDefinition,
26+ opts: { gpu?: boolean } = {},
27+): string {
28+ const g = opts.gpu ? method.gpu : undefined;
29+ const solve = g
30+ ? [
31+ `% Solve. With the GPU option the host runs ${g.setupFunction} and ${g.outputFunction} here in numbl and`,
32+ `% ${g.loopFile} as its WebGPU port; without WebGPU, or when this script or ${g.loopFile} is edited,`,
33+ `% all three run in numbl. The library files are in the other tabs.`,
34+ `S = ${g.setupFunction}(params);`,
35+ `R = ${g.loopFile.replace(/\.m$/, "")}(S); % replaced by the WebGPU port when the loop runs on the GPU`,
36+ `result = ${g.outputFunction}(S, R);`,
37+ ]
38+ : [
39+ `% Solve. ${method.solverFunction} and the functions it calls are the library files in the other tabs.`,
40+ `result = ${method.solverFunction}(params);`,
41+ ];
42+ return [generateParamsBlock(spec), "", ...solve, ""].join("\n");
2643 }
2744
2845 /** Bundle for running one method on a spec. */
29-export function buildRunBundle(spec: ProblemSpec, method: MethodDefinition): SourceBundle {
46+export function buildRunBundle(
47+ spec: ProblemSpec,
48+ method: MethodDefinition,
49+ opts: { gpu?: boolean } = {},
50+): SourceBundle {
3051 return {
31- files: [{ name: MAIN_FILE, source: generateMain(spec, method) }, ...LIBRARY_FILES],
52+ files: [{ name: MAIN_FILE, source: generateMain(spec, method, opts) }, ...LIBRARY_FILES],
3253 entry: MAIN_FILE,
3354 resultVariable: "result",
3455 };
3556 }
3657
58+/**
59+ * The two numbl scripts of a GPU run (see runner/executeGpu.ts): the setup,
60+ * whose struct S the host reads (sparse operators as find() triplets, since
61+ * a sparse matrix does not cross to the host), and the output, which takes the
62+ * loop's record as the workspace variable R. Null for a method without a port.
63+ */
64+export function buildGpuStages(spec: ProblemSpec, method: MethodDefinition): GpuStages | null {
65+ const g = method.gpu;
66+ if (!g) return null;
67+ const params = generateParamsBlock(spec);
68+ const setup = [
69+ params,
70+ "",
71+ `% Stage 1 of a GPU run: the setup, read by the WebGPU port of ${g.loopFile}.`,
72+ `S = ${g.setupFunction}(params);`,
73+ "if S.sparse",
74+ " % the sparse operators cross to the host as triplets",
75+ " [ki, kj, kv] = find(S.Km); S.Km_i = ki; S.Km_j = kj; S.Km_v = kv;",
76+ " [li, lj, lv] = find(S.Lm); S.Lm_i = li; S.Lm_j = lj; S.Lm_v = lv;",
77+ "end",
78+ "",
79+ ].join("\n");
80+ const output = [
81+ params,
82+ "",
83+ `% Stage 3 of a GPU run: R is the record of the WebGPU loop, passed in by the host.`,
84+ `S = ${g.setupFunction}(params);`,
85+ `result = ${g.outputFunction}(S, R);`,
86+ "",
87+ ].join("\n");
88+ return { setup, output, loopFile: g.loopFile };
89+}
90+
3791 /** Bundle that evaluates a method's cost plan (currently only the FD method has one). */
3892 export function buildPlanBundle(spec: ProblemSpec, method: MethodDefinition): SourceBundle | null {
3993 if (!method.planFunction) return null;
40- const main = [generateParamsBlock(spec), "", `plan = ${method.planFunction}(params);`, ""].join("\n");
94+ const main = [generateParamsBlock(spec), "", `plan = ${method.planFunction}(params);`, ""].join(
95+ "\n",
96+ );
4197 return {
4298 files: [{ name: "plan.m", source: main }, ...LIBRARY_FILES],
4399 entry: "plan.m",
@@ -67,6 +123,8 @@ export function buildCompareBundle(): SourceBundle {
67123 export function withEdits(bundle: SourceBundle, edits: Record<string, string>): SourceBundle {
68124 return {
69125 ...bundle,
70- files: bundle.files.map((f) => (edits[f.name] !== undefined ? { ...f, source: edits[f.name] } : f)),
126+ files: bundle.files.map((f) =>
127+ edits[f.name] !== undefined ? { ...f, source: edits[f.name] } : f,
128+ ),
71129 };
72130 }
packages/core/src/generate/library.tsmodified+39−4View file
@@ -1,6 +1,9 @@
11 import plate_modal from "../m/plate_modal.m?raw";
22 import plate_fd from "../m/plate_fd.m?raw";
33 import plate_fd_plan from "../m/plate_fd_plan.m?raw";
4+import plate_fd_setup from "../m/plate_fd_setup.m?raw";
5+import plate_fd_loop from "../m/plate_fd_loop.m?raw";
6+import plate_fd_output from "../m/plate_fd_output.m?raw";
47 import plate_setup from "../m/plate_setup.m?raw";
58 import plate_force_time from "../m/plate_force_time.m?raw";
69 import plate_force_footprint from "../m/plate_force_footprint.m?raw";
@@ -17,6 +20,8 @@ import plate_implicit_plan from "../m/plate_implicit_plan.m?raw";
1720 import plate_matrices from "../m/plate_matrices.m?raw";
1821 import plate_rigid_modes from "../m/plate_rigid_modes.m?raw";
1922 import plate_curvature_free from "../m/plate_curvature_free.m?raw";
23+import plate_domain from "../m/plate_domain.m?raw";
24+import plate_curvature_domain from "../m/plate_curvature_domain.m?raw";
2025 import plate_result_schema from "../m/plate_result_schema.m?raw";
2126
2227 export interface SourceFile {
@@ -36,7 +41,26 @@ export const LIBRARY_FILES: LibraryFile[] = [
3641 source: plate_modal,
3742 summary: "Modal superposition on a simply supported rectangle",
3843 },
39- { name: "plate_fd.m", source: plate_fd, summary: "Explicit finite-difference time stepping" },
44+ {
45+ name: "plate_fd.m",
46+ source: plate_fd,
47+ summary: "Explicit finite differences: setup, time loop and output in turn",
48+ },
49+ {
50+ name: "plate_fd_setup.m",
51+ source: plate_fd_setup,
52+ summary: "Grid, step, load, probe, listener and operators for the explicit loop",
53+ },
54+ {
55+ name: "plate_fd_loop.m",
56+ source: plate_fd_loop,
57+ summary: "Explicit time stepping; the lab has a WebGPU port of this file",
58+ },
59+ {
60+ name: "plate_fd_output.m",
61+ source: plate_fd_output,
62+ summary: "Stresses, audio, spectra and the result struct from the loop's record",
63+ },
4064 {
4165 name: "plate_fd_plan.m",
4266 source: plate_fd_plan,
@@ -52,10 +76,16 @@ export const LIBRARY_FILES: LibraryFile[] = [
5276 source: plate_implicit_plan,
5377 summary: "Grid, accuracy-chosen time step and size of the implicit scheme",
5478 },
79+ {
80+ name: "plate_domain.m",
81+ source: plate_domain,
82+ summary:
83+ "The plate's footprint on the grid: every node for a rectangle, a staircase for an ellipse",
84+ },
5585 {
5686 name: "plate_matrices.m",
5787 source: plate_matrices,
58- summary: "Sparse mass, Laplacian and stiffness matrices for every edge condition",
88+ summary: "Sparse mass, Laplacian and stiffness matrices for every edge condition and shape",
5989 },
6090 {
6191 name: "plate_rigid_modes.m",
@@ -75,7 +105,7 @@ export const LIBRARY_FILES: LibraryFile[] = [
75105 {
76106 name: "plate_force_footprint.m",
77107 source: plate_force_footprint,
78- summary: "Raised-cosine spatial footprint of the load",
108+ summary: "Raised-cosine spatial footprint of the load, restricted to the plate's nodes",
79109 },
80110 {
81111 name: "plate_stress.m",
@@ -107,7 +137,12 @@ export const LIBRARY_FILES: LibraryFile[] = [
107137 {
108138 name: "plate_curvature_free.m",
109139 source: plate_curvature_free,
110- summary: "Curvatures of a free plate, edges from the moment condition",
140+ summary: "Curvatures of a free rectangle, edges from the moment condition",
141+ },
142+ {
143+ name: "plate_curvature_domain.m",
144+ source: plate_curvature_domain,
145+ summary: "Curvatures on an ellipse: central differences inside, copied along the staircase",
111146 },
112147 { name: "plate_compare.m", source: plate_compare, summary: "Error norms between two solutions" },
113148 {
packages/core/src/generate/params.tsmodified+2−2View file
@@ -15,9 +15,9 @@ export function generateParamsBlock(spec: ProblemSpec): string {
1515 "% Problem specification, generated from the controls. SI units: m, s, N, Pa, kg/m^3.",
1616 "% Edit anything here or in the library files; the controls regenerate only this block.",
1717 "params = struct();",
18- `params.plate = struct('Lx', ${n(spec.plate.Lx)}, 'Ly', ${n(spec.plate.Ly)}, 'h', ${n(spec.plate.h)}); % side lengths and thickness [m]`,
18+ `params.plate = struct('shape', ${s(spec.plate.shape)}, 'Lx', ${n(spec.plate.Lx)}, 'Ly', ${n(spec.plate.Ly)}, 'h', ${n(spec.plate.h)}); % 'rectangle' or 'ellipse' in the box Lx x Ly, thickness h [m]`,
1919 `params.material = struct('name', ${s(spec.material.name)}, 'E', ${n(spec.material.E)}, 'nu', ${n(spec.material.nu)}, 'rho', ${n(spec.material.rho)}); % Young's modulus [Pa], Poisson's ratio, density [kg/m^3]`,
20- `params.bc = ${s(spec.bc)}; % 'simply-supported', 'clamped' or 'free' (modal: simply-supported only)`,
20+ `params.bc = ${s(spec.bc)}; % 'simply-supported', 'clamped', 'free' or 'pinned-center' (modal: simply-supported only)`,
2121 `params.damping = struct('sigma0', ${n(spec.damping.sigma0)}, 'sigma1', ${n(spec.damping.sigma1)}); % loss rates: frequency-independent [1/s], frequency-dependent [m^2/s]`,
2222 `params.excitation = ${excitation}`,
2323 `params.probe = struct('x', ${n(spec.probe.x)}, 'y', ${n(spec.probe.y)}); % where the time series are recorded [m]`,
packages/core/src/gpu/fdLoop.tsadded+664−0View file
@@ -0,0 +1,664 @@
1+/**
2+ * WebGPU port of plate_fd_loop.m: the explicit time stepping of the plate,
3+ * run as compute shaders. Everything around the loop stays in numbl: the host
4+ * runs the unmodified plate_fd_setup.m, hands its struct in here as FdLoopSetup,
5+ * and passes the record this returns to the unmodified plate_fd_output.m, which
6+ * computes the stresses, the audio and the spectra in double precision.
7+ *
8+ * One time step is three compute dispatches for a rectangle with supported or
9+ * clamped edges, and three or five when the setup carries sparse operators
10+ * (free edges, or an ellipse of any edge condition); the dispatch boundaries
11+ * are the global synchronization points the stencil needs:
12+ * op Lw = Laplacian of w at every node (ghost layer by reflection), or,
13+ * with sparse operators, Lw = Lm w / mvec and B = Km w / mvec
14+ * update the new field, wn = c_new (2 w - c_old wm - k2dt2 B + s1dt2 (Lw - Lw_m)
15+ * + fdt2 F(n) s); edges zeroed for the dense rectangle, the nodes
16+ * that are not unknowns (S.keep = 0) zeroed for the sparse form
17+ * dots (nr > 0) the rigid-body components Qb' wn
18+ * project (nr > 0) wn - Rb (Qb' wn)
19+ * reduce acceleration summed per listener delay bin into the audio bins,
20+ * the probe displacement and curvatures, max |wn| for the blow-up test
21+ *
22+ * Steps are recorded in chunks of a few thousand into one command buffer each,
23+ * with the step index a chunk base (ctl[0], written by the host once per chunk)
24+ * plus a within-chunk offset from a dynamic-offset uniform; every recorded
25+ * dispatch is identical. Two chunks are in flight at a time and the blow-up
26+ * flag is read back after each, so a diverging run stops a chunk late rather
27+ * than at the end.
28+ *
29+ * Precision. The GPU works in single precision. The state is carried as
30+ * (w, d = w - wm) rather than (w, wm): the update is algebraically the .m's,
31+ * dn = (d + rest) - cs (2 d + rest), wn = w + dn, cs = sigma0 dt / (1 + sigma0 dt),
32+ * with rest = -k2dt2 B + s1dt2 (Lw - Lw_m) + fdt2 F(n) s, so that no factor of
33+ * the form 1 - epsilon is stored in f32 and the rounding of w (about 6e-8 of a
34+ * micrometre-scale displacement) does not feed back into the velocity as it
35+ * would through 2 w - wm. The acceleration for the audio is (dn - d) / dt^2.
36+ */
37+
38+export interface Triplets {
39+ /** 1-based row and column indices, as MATLAB's find returns them */
40+ i: Float64Array;
41+ j: Float64Array;
42+ v: Float64Array;
43+}
44+
45+/** The fields of plate_fd_setup's struct that the loop reads; indices are 1-based as in the .m. */
46+export interface FdLoopSetup {
47+ nx: number;
48+ ny: number;
49+ nsteps: number;
50+ dx: number;
51+ dy: number;
52+ dt: number;
53+ fs: number;
54+ bc: string;
55+ /** Load footprint, ny x nx column-major */
56+ s: Float64Array;
57+ /** Force at the nsteps + 1 time levels */
58+ F: Float64Array;
59+ /** Probe row and column */
60+ ip: number;
61+ jp: number;
62+ /** Steps at which the field is kept, sorted and unique */
63+ frameSteps: Float64Array;
64+ /** Listener delay bin and weight of every node */
65+ bins: Float64Array;
66+ wgt: Float64Array;
67+ nbins: number;
68+ na: number;
69+ /** The loop uses the sparse operators (free edges, or an ellipse) */
70+ sparse: boolean;
71+ /** Sparse form only: 1 at the unknowns and 0 at the nodes held at zero (N), node weights (N) */
72+ keep: Float64Array | null;
73+ mvec: Float64Array | null;
74+ /** Rigid-body motions to remove: basis and projector (N x nr), nr = 0, 2 or 3 */
75+ Rb: Float64Array | null;
76+ Qb: Float64Array | null;
77+ nr: number;
78+ Km: Triplets | null;
79+ Lm: Triplets | null;
80+ sigma0: number;
81+ k2dt2: number;
82+ s1dt2: number;
83+ fdt2: number;
84+ blowup: number;
85+}
86+
87+/** What plate_fd_loop.m returns, in single precision, before conversion for plate_fd_output.m. */
88+export interface FdLoopRecord {
89+ w_p: Float32Array;
90+ wxx_p: Float32Array;
91+ wyy_p: Float32Array;
92+ wxy_p: Float32Array;
93+ /** ny x nx x nf */
94+ Wf: Float32Array;
95+ /** nbins x na */
96+ Pbins: Float32Array;
97+ divergedStep: number;
98+ gpuMs: number;
99+}
100+
101+export interface GpuHandle {
102+ device: GPUDevice;
103+ adapter: string;
104+}
105+
106+export function describeAdapter(adapter: GPUAdapter): string {
107+ const info = (adapter as { info?: Partial<GPUAdapterInfo> }).info;
108+ if (!info) return "unknown adapter";
109+ const name = [info.vendor, info.architecture].filter(Boolean).join(" ");
110+ return name || info.description || "unknown adapter";
111+}
112+
113+/** The device request the browser worker and the Node test both make. */
114+export async function requestGpuDevice(): Promise<GpuHandle> {
115+ const gpu = (globalThis as { navigator?: { gpu?: GPU } }).navigator?.gpu;
116+ if (!gpu) throw new Error("WebGPU is not available");
117+ const adapter = await gpu.requestAdapter();
118+ if (!adapter) throw new Error("no WebGPU adapter");
119+ const device = await adapter.requestDevice({
120+ requiredLimits: {
121+ maxStorageBufferBindingSize: adapter.limits.maxStorageBufferBindingSize,
122+ maxBufferSize: adapter.limits.maxBufferSize,
123+ },
124+ });
125+ return { device, adapter: describeAdapter(adapter) };
126+}
127+
128+const WG = 64;
129+const RWG = 256;
130+const CHUNK = 2048;
131+const NONE = 0xffffffff;
132+
133+const fmt = (x: number): string => {
134+ if (!Number.isFinite(x)) throw new Error("non-finite kernel constant");
135+ let s = x.toPrecision(9);
136+ if (!/[.e]/.test(s)) s += ".0";
137+ return s + "f";
138+};
139+
140+interface Csr {
141+ ptr: Uint32Array;
142+ col: Uint32Array;
143+ val: Float32Array;
144+}
145+
146+function toCsr(t: Triplets, n: number): Csr {
147+ const count = new Uint32Array(n);
148+ for (let p = 0; p < t.i.length; p++) count[t.i[p] - 1]++;
149+ const ptr = new Uint32Array(n + 1);
150+ for (let r = 0; r < n; r++) ptr[r + 1] = ptr[r] + count[r];
151+ const fill = ptr.slice(0, n);
152+ const col = new Uint32Array(t.i.length);
153+ const val = new Float32Array(t.i.length);
154+ for (let p = 0; p < t.i.length; p++) {
155+ const q = fill[t.i[p] - 1]++;
156+ col[q] = t.j[p] - 1;
157+ val[q] = t.v[p];
158+ }
159+ return { ptr, col, val };
160+}
161+
162+/** Nodes grouped by listener delay bin, as a CSR-like list. */
163+function binLists(bins: Float64Array, nbins: number): { ptr: Uint32Array; nodes: Uint32Array } {
164+ const n = bins.length;
165+ const count = new Uint32Array(nbins);
166+ for (let i = 0; i < n; i++) count[bins[i] - 1]++;
167+ const ptr = new Uint32Array(nbins + 1);
168+ for (let b = 0; b < nbins; b++) ptr[b + 1] = ptr[b] + count[b];
169+ const fill = ptr.slice(0, nbins);
170+ const nodes = new Uint32Array(n);
171+ for (let i = 0; i < n; i++) nodes[fill[bins[i] - 1]++] = i;
172+ return { ptr, nodes };
173+}
174+
175+function pack<T extends Uint32Array | Float32Array>(
176+ parts: [string, ArrayLike<number>][],
177+ make: (n: number) => T,
178+): { data: T; off: Record<string, number> } {
179+ const off: Record<string, number> = {};
180+ let cursor = 0;
181+ for (const [name, arr] of parts) {
182+ off[name] = cursor;
183+ cursor += arr.length;
184+ }
185+ const data = make(Math.max(1, cursor));
186+ for (const [name, arr] of parts) data.set(arr, off[name]);
187+ return { data, off };
188+}
189+
190+function buildShaders(su: FdLoopSetup) {
191+ const nx = su.nx | 0;
192+ const ny = su.ny | 0;
193+ const N = nx * ny;
194+ const nsteps = su.nsteps | 0;
195+ const nbins = su.nbins | 0;
196+ const na = su.na | 0;
197+ const nf = su.frameSteps.length;
198+ const sp = su.sparse;
199+ const nr = sp ? su.nr : 0;
200+ // With a projection the frames are written after it; otherwise straight from the update
201+ const framesInUpdate = nr === 0;
202+ if (su.s.length !== N || su.bins.length !== N || su.wgt.length !== N)
203+ throw new Error("setup arrays do not match the grid");
204+ if (su.F.length < nsteps) throw new Error("force history shorter than the step count");
205+
206+ // ---- per-step tables, decided in f64 exactly as the .m decides them
207+ const frameOf = new Uint32Array(nsteps).fill(NONE);
208+ for (let kf = 0; kf < nf; kf++) {
209+ const step = su.frameSteps[kf];
210+ if (step >= 1 && step <= nsteps) frameOf[step - 1] = kf;
211+ }
212+ const iaOf = new Uint32Array(nsteps);
213+ for (let k = 0; k < nsteps; k++) {
214+ const nstep = k + 1;
215+ const ia = Math.floor(nstep * su.dt * su.fs); // ia = floor(nstep * dt * fs) + 1, made 0-based
216+ if (ia >= na) throw new Error("audio bin index out of range");
217+ iaOf[k] = ia;
218+ }
219+ const bl = binLists(su.bins, nbins);
220+
221+ const f32 = (a: ArrayLike<number>) => Float32Array.from(a);
222+ const opsParts: [string, ArrayLike<number>][] = [
223+ ["S", f32(su.s)],
224+ ["WGT", f32(su.wgt)],
225+ ["F", f32(su.F.subarray(0, nsteps))],
226+ ];
227+ const idxParts: [string, ArrayLike<number>][] = [
228+ ["BINPTR", bl.ptr],
229+ ["BINNODES", bl.nodes],
230+ ["FRAME", frameOf],
231+ ["IA", iaOf],
232+ ];
233+ if (sp) {
234+ if (!su.mvec || !su.keep || !su.Km || !su.Lm)
235+ throw new Error("the sparse form needs mvec, keep, Km and Lm");
236+ if (su.mvec.length !== N || su.keep.length !== N)
237+ throw new Error("node arrays do not match the grid");
238+ const K = toCsr(su.Km, N);
239+ const L = toCsr(su.Lm, N);
240+ opsParts.push(["MVEC", f32(su.mvec)], ["KEEP", f32(su.keep)], ["KV", K.val], ["LV", L.val]);
241+ idxParts.push(["KPTR", K.ptr], ["KCOL", K.col], ["LPTR", L.ptr], ["LCOL", L.col]);
242+ if (nr > 0) {
243+ if (!su.Rb || !su.Qb || su.Rb.length !== N * nr || su.Qb.length !== N * nr)
244+ throw new Error("rigid-body arrays do not match the grid");
245+ opsParts.push(["RB", f32(su.Rb)], ["QB", f32(su.Qb)]);
246+ }
247+ }
248+ const ops = pack(opsParts, (n) => new Float32Array(n));
249+ const idx = pack(idxParts, (n) => new Uint32Array(n));
250+
251+ // ---- state buffer: two slots each of w, d = w - wm and Lw; B; a few scalars
252+ const ST_W = 0;
253+ const ST_D = 2 * N;
254+ const ST_L = 4 * N;
255+ const ST_B = 6 * N;
256+ const ST_AUX = 7 * N;
257+ const stLen = ST_AUX + 16;
258+
259+ const nprobe = nsteps + 1;
260+ const OUT_PB = 4 * nprobe;
261+ const outLen = OUT_PB + nbins * na;
262+
263+ const s0dt = su.sigma0 * su.dt;
264+ const c: Record<string, number> = {
265+ INV_DX2: 1 / (su.dx * su.dx),
266+ INV_DY2: 1 / (su.dy * su.dy),
267+ INV_DXDY4: 1 / (4 * su.dx * su.dy),
268+ K2DT2: su.k2dt2,
269+ S1DT2: su.s1dt2,
270+ FDT2: su.fdt2,
271+ CS: s0dt / (1 + s0dt),
272+ INV_DT2: 1 / (su.dt * su.dt),
273+ DTFS: su.dt * su.fs,
274+ BLOWUP: su.blowup,
275+ SGN: su.bc === "clamped" ? 1 : -1,
276+ };
277+ const ip = su.ip - 1;
278+ const jp = su.jp - 1;
279+ if (ip < 1 || ip > ny - 2 || jp < 1 || jp > nx - 2)
280+ throw new Error("probe is not an interior node");
281+
282+ const header = `
283+const NX = ${nx}u; const NY = ${ny}u; const N = ${N}u;
284+const NXI: i32 = ${nx}; const NYI: i32 = ${ny};
285+const NSTEPS = ${nsteps}u; const NPROBE = ${nprobe}u;
286+const NBINS = ${nbins}u; const NA = ${na}u;
287+const IPROBE = ${ip + jp * ny}u;
288+const NR = ${nr}u;
289+${Object.entries(c)
290+ .map(([k, v]) => `const ${k}: f32 = ${fmt(v)};`)
291+ .join("\n")}
292+${Object.entries(ops.off)
293+ .map(([k, v]) => `const OP_${k} = ${v}u;`)
294+ .join("\n")}
295+${Object.entries(idx.off)
296+ .map(([k, v]) => `const IX_${k} = ${v}u;`)
297+ .join("\n")}
298+const ST_W = ${ST_W}u; const ST_D = ${ST_D}u; const ST_L = ${ST_L}u; const ST_B = ${ST_B}u; const ST_AUX = ${ST_AUX}u;
299+const OUT_PB = ${OUT_PB}u;
300+
301+@group(0) @binding(0) var<storage, read> ops: array<f32>;
302+@group(0) @binding(1) var<storage, read> idx: array<u32>;
303+@group(0) @binding(2) var<storage, read_write> st: array<f32>;
304+@group(0) @binding(3) var<storage, read_write> frames: array<f32>;
305+@group(0) @binding(4) var<storage, read_write> outb: array<f32>;
306+@group(0) @binding(5) var<storage, read_write> ctl: array<u32>;
307+@group(0) @binding(6) var<uniform> uni: vec4<u32>;
308+
309+// The blow-up flag holds the 1-based step at which |w| passed the limit; the
310+// steps after it do nothing, the step itself still records, as the .m does.
311+fn stopped(k: u32) -> bool {
312+ let d = ctl[1];
313+ return d != 0u && d <= k;
314+}
315+// Slot of the current field (w, d, Lw) at step k and of the one being written
316+fn cur(k: u32) -> u32 { return (k & 1u) * N; }
317+fn nxt(k: u32) -> u32 { return ((k + 1u) & 1u) * N; }
318+`;
319+
320+ // Padded field of plate_pad: the edge nodes are zero and one ghost layer
321+ // reflects the first interior value, odd for supported, even for clamped.
322+ const opSupported = `${header}
323+fn wat(wS: u32, i: i32, j: i32) -> f32 {
324+ var ii = i; var jj = j; var sg: f32 = 1.0;
325+ if (ii < 0) { ii = 1; sg = sg * SGN; } else if (ii >= NYI) { ii = NYI - 2; sg = sg * SGN; }
326+ if (jj < 0) { jj = 1; sg = sg * SGN; } else if (jj >= NXI) { jj = NXI - 2; sg = sg * SGN; }
327+ return sg * st[wS + u32(ii) + u32(jj) * NY];
328+}
329+@compute @workgroup_size(${WG})
330+fn main(@builtin(global_invocation_id) gid: vec3<u32>) {
331+ let k = ctl[0] + uni.x;
332+ if (k >= NSTEPS) { return; }
333+ let i = gid.x;
334+ if (i >= N) { return; }
335+ if (stopped(k)) { return; }
336+ let wS = ST_W + cur(k);
337+ let r = i32(i % NY); let cc = i32(i / NY);
338+ let w0 = st[wS + i];
339+ let lw = (wat(wS, r + 1, cc) - 2.0 * w0 + wat(wS, r - 1, cc)) * INV_DY2
340+ + (wat(wS, r, cc + 1) - 2.0 * w0 + wat(wS, r, cc - 1)) * INV_DX2;
341+ st[ST_L + cur(k) + i] = lw;
342+}
343+`;
344+
345+ // Sparse form: Lw = Lm w / mvec and B = Km w / mvec from the sparse operators
346+ const opSparse = `${header}
347+@compute @workgroup_size(${WG})
348+fn main(@builtin(global_invocation_id) gid: vec3<u32>) {
349+ let k = ctl[0] + uni.x;
350+ if (k >= NSTEPS) { return; }
351+ let i = gid.x;
352+ if (i >= N) { return; }
353+ if (stopped(k)) { return; }
354+ let wS = ST_W + cur(k);
355+ var lsum: f32 = 0.0;
356+ for (var p = idx[IX_LPTR + i]; p < idx[IX_LPTR + i + 1u]; p = p + 1u) {
357+ lsum = lsum + ops[OP_LV + p] * st[wS + idx[IX_LCOL + p]];
358+ }
359+ var ksum: f32 = 0.0;
360+ for (var p = idx[IX_KPTR + i]; p < idx[IX_KPTR + i + 1u]; p = p + 1u) {
361+ ksum = ksum + ops[OP_KV + p] * st[wS + idx[IX_KCOL + p]];
362+ }
363+ let m = ops[OP_MVEC + i];
364+ st[ST_L + cur(k) + i] = lsum / m;
365+ st[ST_B + i] = ksum / m;
366+}
367+`;
368+
369+ const biharmonic = sp
370+ ? `let b = st[ST_B + i];`
371+ : `let lS = ST_L + cur(k);
372+ let lw0 = st[lS + i];
373+ let b = (st[lS + i + 1u] - 2.0 * lw0 + st[lS + i - 1u]) * INV_DY2
374+ + (st[lS + i + NY] - 2.0 * lw0 + st[lS + i - NY]) * INV_DX2;`;
375+ const update = `${header}
376+@compute @workgroup_size(${WG})
377+fn main(@builtin(global_invocation_id) gid: vec3<u32>) {
378+ let k = ctl[0] + uni.x;
379+ if (k >= NSTEPS) { return; }
380+ let i = gid.x;
381+ if (i >= N) { return; }
382+ if (stopped(k)) { return; }
383+ let r = i % NY; let cc = i / NY;
384+ let wS = ST_W + cur(k); let wnS = ST_W + nxt(k);
385+ let dS = ST_D + cur(k); let dnS = ST_D + nxt(k);
386+ var wn: f32 = 0.0; var dn: f32 = 0.0;
387+ ${sp ? "if (true) {" : "if (r > 0u && r < NY - 1u && cc > 0u && cc < NX - 1u) {"}
388+ ${biharmonic}
389+ let lw = st[ST_L + cur(k) + i];
390+ let lwm = st[ST_L + nxt(k) + i];
391+ let rest = -K2DT2 * b + S1DT2 * (lw - lwm) + FDT2 * ops[OP_F + k] * ops[OP_S + i];
392+ let d = st[dS + i];
393+ dn = (d + rest) - CS * (2.0 * d + rest);
394+ wn = st[wS + i] + dn;
395+ }
396+ ${sp ? "let kp = ops[OP_KEEP + i]; wn = wn * kp; dn = dn * kp;" : ""}
397+ st[wnS + i] = wn;
398+ st[dnS + i] = dn;
399+ ${
400+ framesInUpdate
401+ ? `let f = idx[IX_FRAME + k];
402+ if (f != ${NONE}u) { frames[f * N + i] = wn; }`
403+ : ""
404+ }
405+}
406+`;
407+
408+ // Rigid-body components Qb' wn (the nodes held at zero already are)
409+ const dots = `${header}
410+var<workgroup> sh0: array<f32, ${RWG}>;
411+var<workgroup> sh1: array<f32, ${RWG}>;
412+var<workgroup> sh2: array<f32, ${RWG}>;
413+@compute @workgroup_size(${RWG})
414+fn main(@builtin(local_invocation_index) lid: u32) {
415+ let k = ctl[0] + uni.x;
416+ let wnS = ST_W + nxt(k);
417+ var a0: f32 = 0.0; var a1: f32 = 0.0; var a2: f32 = 0.0;
418+ for (var i = lid; i < N; i = i + ${RWG}u) {
419+ let wn = st[wnS + i];
420+ a0 = a0 + ops[OP_QB + i] * wn;
421+ a1 = a1 + ops[OP_QB + N + i] * wn;
422+ if (NR == 3u) { a2 = a2 + ops[OP_QB + 2u * N + i] * wn; }
423+ }
424+ sh0[lid] = a0; sh1[lid] = a1; sh2[lid] = a2;
425+ workgroupBarrier();
426+ for (var s = ${RWG / 2}u; s > 0u; s = s >> 1u) {
427+ if (lid < s) { sh0[lid] = sh0[lid] + sh0[lid + s]; sh1[lid] = sh1[lid] + sh1[lid + s]; sh2[lid] = sh2[lid] + sh2[lid + s]; }
428+ workgroupBarrier();
429+ }
430+ if (lid == 0u) { st[ST_AUX] = sh0[0]; st[ST_AUX + 1u] = sh1[0]; st[ST_AUX + 2u] = sh2[0]; }
431+}
432+`;
433+
434+ const project = `${header}
435+@compute @workgroup_size(${WG})
436+fn main(@builtin(global_invocation_id) gid: vec3<u32>) {
437+ let k = ctl[0] + uni.x;
438+ if (k >= NSTEPS) { return; }
439+ let i = gid.x;
440+ if (i >= N) { return; }
441+ if (stopped(k)) { return; }
442+ let wnS = ST_W + nxt(k); let dnS = ST_D + nxt(k);
443+ var corr = ops[OP_RB + i] * st[ST_AUX] + ops[OP_RB + N + i] * st[ST_AUX + 1u];
444+ if (NR == 3u) { corr = corr + ops[OP_RB + 2u * N + i] * st[ST_AUX + 2u]; }
445+ // Rb is zero at the nodes held at zero; the factor makes that exact in f32 as well
446+ let kp = ops[OP_KEEP + i];
447+ let wn = (st[wnS + i] - corr) * kp; let dn = (st[dnS + i] - corr) * kp;
448+ st[wnS + i] = wn; st[dnS + i] = dn;
449+ let f = idx[IX_FRAME + k];
450+ if (f != ${NONE}u) { frames[f * N + i] = wn; }
451+}
452+`;
453+
454+ // Workgroups 0..NBINS-1 sum the weighted acceleration of their bin's nodes into
455+ // the audio sample of this step; workgroup NBINS takes max |wn| for the
456+ // blow-up test and records the probe.
457+ const reduce = `${header}
458+var<workgroup> sh: array<f32, ${RWG}>;
459+@compute @workgroup_size(${RWG})
460+fn main(@builtin(local_invocation_index) lid: u32, @builtin(workgroup_id) wgid: vec3<u32>) {
461+ let k = ctl[0] + uni.x;
462+ let stop = stopped(k) || k >= NSTEPS;
463+ let wnS = ST_W + nxt(k);
464+ let dS = ST_D + cur(k); let dnS = ST_D + nxt(k);
465+ let wg = wgid.x;
466+ if (wg < NBINS) {
467+ let p0 = idx[IX_BINPTR + wg]; let p1 = idx[IX_BINPTR + wg + 1u];
468+ var a: f32 = 0.0;
469+ for (var p = p0 + lid; p < p1; p = p + ${RWG}u) {
470+ let i = idx[IX_BINNODES + p];
471+ a = a + ops[OP_WGT + i] * ((st[dnS + i] - st[dS + i]) * INV_DT2);
472+ }
473+ sh[lid] = a;
474+ workgroupBarrier();
475+ for (var s = ${RWG / 2}u; s > 0u; s = s >> 1u) {
476+ if (lid < s) { sh[lid] = sh[lid] + sh[lid + s]; }
477+ workgroupBarrier();
478+ }
479+ if (lid == 0u && !stop) {
480+ let ia = idx[IX_IA + k];
481+ let o = OUT_PB + wg + ia * NBINS;
482+ outb[o] = outb[o] + sh[0] * DTFS;
483+ }
484+ } else {
485+ var m: f32 = 0.0;
486+ for (var i = lid; i < N; i = i + ${RWG}u) { m = max(m, abs(st[wnS + i])); }
487+ sh[lid] = m;
488+ workgroupBarrier();
489+ for (var s = ${RWG / 2}u; s > 0u; s = s >> 1u) {
490+ if (lid < s) { sh[lid] = max(sh[lid], sh[lid + s]); }
491+ workgroupBarrier();
492+ }
493+ if (lid == 0u && !stop) {
494+ let c = st[wnS + IPROBE];
495+ let e = st[wnS + IPROBE + NY]; let w = st[wnS + IPROBE - NY];
496+ let n = st[wnS + IPROBE + 1u]; let so = st[wnS + IPROBE - 1u];
497+ let ne = st[wnS + IPROBE + NY + 1u]; let nw = st[wnS + IPROBE - NY + 1u];
498+ let se = st[wnS + IPROBE + NY - 1u]; let sw = st[wnS + IPROBE - NY - 1u];
499+ outb[k + 1u] = c;
500+ outb[NPROBE + k + 1u] = (e - 2.0 * c + w) * INV_DX2;
501+ outb[2u * NPROBE + k + 1u] = (n - 2.0 * c + so) * INV_DY2;
502+ outb[3u * NPROBE + k + 1u] = (ne - se - nw + sw) * INV_DXDY4;
503+ if (sh[0] >= BLOWUP) { ctl[1] = k + 1u; }
504+ }
505+ }
506+}
507+`;
508+
509+ const op = sp ? opSparse : opSupported;
510+ const sources = nr > 0 ? [op, update, dots, project, reduce] : [op, update, reduce];
511+ const groups =
512+ nr > 0
513+ ? [Math.ceil(N / WG), Math.ceil(N / WG), 1, Math.ceil(N / WG), nbins + 1]
514+ : [Math.ceil(N / WG), Math.ceil(N / WG), nbins + 1];
515+ return { N, nf, nprobe, nbins, na, outLen, stLen, ops: ops.data, idx: idx.data, sources, groups };
516+}
517+
518+export async function runFdLoopGpu(
519+ device: GPUDevice,
520+ su: FdLoopSetup,
521+ onProgress?: (stepsDone: number, nsteps: number) => void,
522+): Promise<FdLoopRecord> {
523+ const nsteps = su.nsteps | 0;
524+ const b = buildShaders(su);
525+
526+ const bgl = device.createBindGroupLayout({
527+ entries: [
528+ { binding: 0, visibility: GPUShaderStage.COMPUTE, buffer: { type: "read-only-storage" } },
529+ { binding: 1, visibility: GPUShaderStage.COMPUTE, buffer: { type: "read-only-storage" } },
530+ { binding: 2, visibility: GPUShaderStage.COMPUTE, buffer: { type: "storage" } },
531+ { binding: 3, visibility: GPUShaderStage.COMPUTE, buffer: { type: "storage" } },
532+ { binding: 4, visibility: GPUShaderStage.COMPUTE, buffer: { type: "storage" } },
533+ { binding: 5, visibility: GPUShaderStage.COMPUTE, buffer: { type: "storage" } },
534+ {
535+ binding: 6,
536+ visibility: GPUShaderStage.COMPUTE,
537+ buffer: { type: "uniform", hasDynamicOffset: true, minBindingSize: 16 },
538+ },
539+ ],
540+ });
541+ const layout = device.createPipelineLayout({ bindGroupLayouts: [bgl] });
542+ const pipelines: GPUComputePipeline[] = [];
543+ for (const code of b.sources) {
544+ const module = device.createShaderModule({ code });
545+ const info = await module.getCompilationInfo();
546+ const errs = info.messages.filter((m) => m.type === "error");
547+ if (errs.length)
548+ throw new Error(
549+ "WGSL: " + errs.map((m) => `${m.lineNum}:${m.linePos} ${m.message}`).join("; "),
550+ );
551+ pipelines.push(
552+ device.createComputePipeline({ layout, compute: { module, entryPoint: "main" } }),
553+ );
554+ }
555+
556+ const mk = (size: number, usage: number) =>
557+ device.createBuffer({ size: Math.max(16, size), usage });
558+ const opsBuf = mk(b.ops.byteLength, GPUBufferUsage.STORAGE | GPUBufferUsage.COPY_DST);
559+ device.queue.writeBuffer(opsBuf, 0, b.ops);
560+ const idxBuf = mk(b.idx.byteLength, GPUBufferUsage.STORAGE | GPUBufferUsage.COPY_DST);
561+ device.queue.writeBuffer(idxBuf, 0, b.idx);
562+ const stBuf = mk(b.stLen * 4, GPUBufferUsage.STORAGE | GPUBufferUsage.COPY_DST);
563+ const framesBuf = mk(b.N * b.nf * 4, GPUBufferUsage.STORAGE | GPUBufferUsage.COPY_SRC);
564+ const outBuf = mk(b.outLen * 4, GPUBufferUsage.STORAGE | GPUBufferUsage.COPY_SRC);
565+ const ctlBuf = mk(16, GPUBufferUsage.STORAGE | GPUBufferUsage.COPY_SRC | GPUBufferUsage.COPY_DST);
566+ const uniBuf = mk(CHUNK * 256, GPUBufferUsage.UNIFORM | GPUBufferUsage.COPY_DST);
567+ {
568+ const uni = new Uint32Array(CHUNK * 64);
569+ for (let s = 0; s < CHUNK; s++) uni[s * 64] = s;
570+ device.queue.writeBuffer(uniBuf, 0, uni);
571+ }
572+ const ctlStage = [0, 1].map(() => mk(16, GPUBufferUsage.MAP_READ | GPUBufferUsage.COPY_DST));
573+ const readFrames = mk(b.N * b.nf * 4, GPUBufferUsage.MAP_READ | GPUBufferUsage.COPY_DST);
574+ const readOut = mk(b.outLen * 4, GPUBufferUsage.MAP_READ | GPUBufferUsage.COPY_DST);
575+ const all = [
576+ opsBuf,
577+ idxBuf,
578+ stBuf,
579+ framesBuf,
580+ outBuf,
581+ ctlBuf,
582+ uniBuf,
583+ ...ctlStage,
584+ readFrames,
585+ readOut,
586+ ];
587+
588+ const bindGroup = device.createBindGroup({
589+ layout: bgl,
590+ entries: [
591+ { binding: 0, resource: { buffer: opsBuf } },
592+ { binding: 1, resource: { buffer: idxBuf } },
593+ { binding: 2, resource: { buffer: stBuf } },
594+ { binding: 3, resource: { buffer: framesBuf } },
595+ { binding: 4, resource: { buffer: outBuf } },
596+ { binding: 5, resource: { buffer: ctlBuf } },
597+ { binding: 6, resource: { buffer: uniBuf, offset: 0, size: 16 } },
598+ ],
599+ });
600+
601+ try {
602+ const t0 = performance.now();
603+ let divergedStep = 0;
604+ const pending: { done: Promise<undefined>; stage: GPUBuffer; kend: number }[] = [];
605+ const settle = async () => {
606+ const p = pending.shift()!;
607+ await p.done;
608+ await p.stage.mapAsync(GPUMapMode.READ);
609+ const flag = new Uint32Array(p.stage.getMappedRange().slice(0))[1];
610+ p.stage.unmap();
611+ onProgress?.(p.kend, nsteps);
612+ if (flag !== 0) divergedStep = flag;
613+ };
614+ let kbase = 0;
615+ let c = 0;
616+ while (kbase < nsteps && divergedStep === 0) {
617+ const n = Math.min(CHUNK, nsteps - kbase);
618+ device.queue.writeBuffer(ctlBuf, 0, Uint32Array.of(kbase));
619+ const enc = device.createCommandEncoder();
620+ const pass = enc.beginComputePass();
621+ for (let s = 0; s < n; s++) {
622+ pass.setBindGroup(0, bindGroup, [s * 256]);
623+ for (let q = 0; q < pipelines.length; q++) {
624+ pass.setPipeline(pipelines[q]);
625+ pass.dispatchWorkgroups(b.groups[q]);
626+ }
627+ }
628+ pass.end();
629+ const stage = ctlStage[c % 2];
630+ enc.copyBufferToBuffer(ctlBuf, 0, stage, 0, 16);
631+ device.queue.submit([enc.finish()]);
632+ pending.push({ done: device.queue.onSubmittedWorkDone(), stage, kend: kbase + n });
633+ kbase += n;
634+ c++;
635+ if (pending.length >= 2) await settle();
636+ }
637+ while (pending.length) await settle();
638+
639+ const enc = device.createCommandEncoder();
640+ enc.copyBufferToBuffer(framesBuf, 0, readFrames, 0, Math.max(16, b.N * b.nf * 4));
641+ enc.copyBufferToBuffer(outBuf, 0, readOut, 0, Math.max(16, b.outLen * 4));
642+ device.queue.submit([enc.finish()]);
643+ await readFrames.mapAsync(GPUMapMode.READ);
644+ await readOut.mapAsync(GPUMapMode.READ);
645+ const Wf = new Float32Array(readFrames.getMappedRange().slice(0)).subarray(0, b.N * b.nf);
646+ const out = new Float32Array(readOut.getMappedRange().slice(0));
647+ readFrames.unmap();
648+ readOut.unmap();
649+ const gpuMs = performance.now() - t0;
650+ const np = b.nprobe;
651+ return {
652+ w_p: out.slice(0, np),
653+ wxx_p: out.slice(np, 2 * np),
654+ wyy_p: out.slice(2 * np, 3 * np),
655+ wxy_p: out.slice(3 * np, 4 * np),
656+ Wf: Float32Array.from(Wf),
657+ Pbins: out.slice(4 * np, 4 * np + b.nbins * b.na),
658+ divergedStep,
659+ gpuMs,
660+ };
661+ } finally {
662+ for (const buf of all) buf.destroy();
663+ }
664+}
packages/core/src/gpu/fdStage.tsadded+126−0View file
@@ -0,0 +1,126 @@
1+/**
2+ * The two hand-offs of a GPU run: plate_fd_setup's struct (as numbl returns it,
3+ * with the sparse operators already turned into find() triplets by the stage
4+ * script) into the loop's input, and the loop's record into the struct R that
5+ * plate_fd_output.m takes. Both sides keep the .m's names and 1-based indices.
6+ */
7+import { isPlainStruct, isPlainTensor, type PlainStruct, type PlainValue } from "../runner/plain";
8+import type { FdLoopRecord, FdLoopSetup, Triplets } from "./fdLoop";
9+
10+class SetupShapeError extends Error {}
11+
12+function field(s: PlainStruct, name: string): PlainValue {
13+ const v = s.fields[name];
14+ if (v === undefined) throw new SetupShapeError(`setup struct is missing "${name}"`);
15+ return v;
16+}
17+
18+function num(s: PlainStruct, name: string): number {
19+ const v = field(s, name);
20+ if (typeof v === "number") return v;
21+ if (typeof v === "boolean") return v ? 1 : 0;
22+ if (isPlainTensor(v) && v.data.length === 1) return v.data[0];
23+ throw new SetupShapeError(`"${name}" should be a scalar`);
24+}
25+
26+function vec(s: PlainStruct, name: string): Float64Array {
27+ const v = field(s, name);
28+ if (typeof v === "number") return Float64Array.of(v);
29+ if (isPlainTensor(v)) return v.data;
30+ throw new SetupShapeError(`"${name}" should be numeric`);
31+}
32+
33+function str(s: PlainStruct, name: string): string {
34+ const v = field(s, name);
35+ if (typeof v === "string") return v;
36+ throw new SetupShapeError(`"${name}" should be a string`);
37+}
38+
39+function sub(s: PlainStruct, name: string): PlainStruct {
40+ const v = field(s, name);
41+ if (!isPlainStruct(v)) throw new SetupShapeError(`"${name}" should be a struct`);
42+ return v;
43+}
44+
45+function triplets(s: PlainStruct, name: string): Triplets {
46+ return { i: vec(s, `${name}_i`), j: vec(s, `${name}_j`), v: vec(s, `${name}_v`) };
47+}
48+
49+/** Read the loop's input from the struct the GPU setup script leaves in S. */
50+export function fdLoopSetupFromPlain(v: PlainValue): FdLoopSetup {
51+ if (!isPlainStruct(v)) throw new SetupShapeError("the setup variable is not a struct");
52+ const P = sub(v, "P");
53+ const sparse = num(v, "sparse") !== 0;
54+ const nr = sparse ? num(v, "nr") : 0;
55+ const nx = num(v, "nx");
56+ const ny = num(v, "ny");
57+ return {
58+ nx,
59+ ny,
60+ nsteps: num(v, "nsteps"),
61+ dx: num(v, "dx"),
62+ dy: num(v, "dy"),
63+ dt: num(v, "dt"),
64+ fs: num(v, "fs"),
65+ bc: str(P, "bc"),
66+ s: vec(v, "s"),
67+ F: vec(v, "F"),
68+ ip: num(v, "ip"),
69+ jp: num(v, "jp"),
70+ frameSteps: vec(v, "frame_steps"),
71+ bins: vec(v, "bins"),
72+ wgt: vec(v, "wgt"),
73+ nbins: num(sub(v, "Lis"), "nbins"),
74+ na: num(v, "na"),
75+ sparse,
76+ keep: sparse ? vec(v, "keep") : null,
77+ mvec: sparse ? vec(v, "mvec") : null,
78+ Rb: nr > 0 ? vec(v, "Rb") : null,
79+ Qb: nr > 0 ? vec(v, "Qb") : null,
80+ nr,
81+ Km: sparse ? triplets(v, "Km") : null,
82+ Lm: sparse ? triplets(v, "Lm") : null,
83+ sigma0: num(sub(P, "damping"), "sigma0"),
84+ k2dt2: num(v, "k2dt2"),
85+ s1dt2: num(v, "s1dt2"),
86+ fdt2: num(v, "fdt2"),
87+ blowup: num(v, "blowup"),
88+ };
89+}
90+
91+/** The record as plate_fd_loop.m would return it, for plate_fd_output.m. */
92+export function fdLoopRecordToPlain(
93+ rec: FdLoopRecord,
94+ su: FdLoopSetup,
95+ adapter: string,
96+): PlainStruct {
97+ const row = (a: Float32Array): PlainValue => ({
98+ kind: "tensor",
99+ shape: [1, a.length],
100+ data: Float64Array.from(a),
101+ });
102+ const nf = su.frameSteps.length;
103+ const notes: string[] = [];
104+ if (rec.divergedStep > 0)
105+ notes.push(
106+ `solution diverged at step ${rec.divergedStep} (t = ${(rec.divergedStep * su.dt).toPrecision(3)} s); results are truncated there`,
107+ );
108+ notes.push(
109+ `time loop ran on the GPU in single precision (${adapter}): ${su.nsteps.toLocaleString()} steps in ${(rec.gpuMs / 1000).toFixed(2)} s`,
110+ );
111+ return {
112+ kind: "struct",
113+ fields: {
114+ w_p: row(rec.w_p),
115+ wxx_p: row(rec.wxx_p),
116+ wyy_p: row(rec.wyy_p),
117+ wxy_p: row(rec.wxy_p),
118+ Wf: { kind: "tensor", shape: [su.ny, su.nx, nf], data: Float64Array.from(rec.Wf) },
119+ Pbins: { kind: "tensor", shape: [su.nbins, su.na], data: Float64Array.from(rec.Pbins) },
120+ diverged_step: rec.divergedStep,
121+ notes: { kind: "cell", shape: [1, notes.length], items: notes },
122+ backend: "webgpu",
123+ loop_s: rec.gpuMs / 1000,
124+ },
125+ };
126+}
packages/core/src/index.tsmodified+8−1View file
@@ -7,7 +7,12 @@ export * from "./methods/registry";
77 export { MODAL_METHOD } from "./methods/modal";
88 export { FD_METHOD } from "./methods/fd";
99 export { IMPLICIT_METHOD } from "./methods/implicit";
10-export { GOVERNING_EQUATIONS, STRESS_EQUATIONS, RADIATION_EQUATIONS } from "./methods/common";
10+export {
11+ GOVERNING_EQUATIONS,
12+ STRESS_EQUATIONS,
13+ RADIATION_EQUATIONS,
14+ CURVED_EDGE_EQUATIONS,
15+} from "./methods/common";
1116 export * from "./generate/format";
1217 export * from "./generate/params";
1318 export * from "./generate/library";
@@ -15,10 +20,12 @@ export * from "./generate/bundle";
1520 export * from "./runner/plain";
1621 export * from "./runner/types";
1722 export { WorkerRunner } from "./runner/workerRunner";
23+export { describeAdapter, type GpuHandle } from "./gpu/fdLoop";
1824 export * from "./results/types";
1925 export * from "./results/parse";
2026 export * from "./fields/field";
2127 export * from "./fields/colormap";
2228 export * from "./fields/raster";
2329 export * from "./audio/buffer";
30+export * from "./audio/wav";
2431 export * from "./cost";
packages/core/src/m/plate_compare.mmodified+2−1View file
@@ -52,7 +52,8 @@ if same_grid
5252 cmp.frames_rel_l2 = zeros(1, nf);
5353 for k = 1:nf
5454 da = A.frames.w(:, :, k); db = B.frames.w(:, :, k);
55- cmp.frames_rel_l2(k) = norm(da(:) - db(:)) / max(norm(da(:)), eps);
55+ fin = isfinite(da) & isfinite(db); % nodes off an elliptical plate are NaN
56+ cmp.frames_rel_l2(k) = norm(da(fin) - db(fin)) / max(norm(da(fin)), eps);
5657 end
5758 cmp.frames_diff = A.frames.w - B.frames.w;
5859 else
packages/core/src/m/plate_curvature_domain.madded+36−0View file
@@ -0,0 +1,36 @@
1+function [wxx, wyy, wxy] = plate_curvature_domain(w, dom)
2+% PLATE_CURVATURE_DOMAIN Second derivatives of the displacement on an embedded plate.
3+%
4+% [wxx, wyy, wxy] = plate_curvature_domain(w, dom) returns the curvatures of
5+% the ny x nx field w on the footprint dom from plate_domain (rows index y,
6+% columns index x). At every node whose eight neighbours belong to the plate
7+% (dom.full3) they are the usual central differences. The remaining nodes of
8+% the plate form the staircase ring along a curved edge, where no single
9+% formula is reliable: an edge node of the staircase has neighbours outside
10+% the plate in directions that depend on where the true edge happens to cut
11+% the grid. Those nodes take the curvatures of the nearest node with a full
12+% stencil (dom.near), a first-order extrapolation that matches the accuracy
13+% of the staircase itself. Outside the plate the curvatures are NaN, so the
14+% stress fields drawn from them are blank there.
15+%
16+% For a clamped edge the largest stress sits on the edge, and this puts it
17+% about one grid spacing inside; refine the grid to sharpen it. The
18+% rectangle keeps its own edge formulas (plate_curvature, plate_curvature_free).
19+
20+dx = dom.dx; dy = dom.dy;
21+[ny, nx] = size(w);
22+wxx = zeros(ny, nx); wyy = zeros(ny, nx); wxy = zeros(ny, nx);
23+c = w(2:end - 1, 2:end - 1);
24+wxx(2:end - 1, 2:end - 1) = (w(2:end - 1, 3:end) - 2 * c + w(2:end - 1, 1:end - 2)) / dx^2;
25+wyy(2:end - 1, 2:end - 1) = (w(3:end, 2:end - 1) - 2 * c + w(1:end - 2, 2:end - 1)) / dy^2;
26+wxy(2:end - 1, 2:end - 1) = (w(3:end, 3:end) - w(3:end, 1:end - 2) - w(1:end - 2, 3:end) + w(1:end - 2, 1:end - 2)) / (4 * dx * dy);
27+ring = find(dom.inside & ~dom.full3);
28+src = dom.near(ring);
29+wxx(ring) = wxx(src);
30+wyy(ring) = wyy(src);
31+wxy(ring) = wxy(src);
32+out = ~dom.inside;
33+wxx(out) = NaN;
34+wyy(out) = NaN;
35+wxy(out) = NaN;
36+end
packages/core/src/m/plate_domain.madded+135−0View file
@@ -0,0 +1,135 @@
1+function dom = plate_domain(plate, nx, ny)
2+% PLATE_DOMAIN The footprint of the plate on a uniform nx x ny grid.
3+%
4+% dom = plate_domain(plate, nx, ny) lays the grid over the bounding box
5+% [0, Lx] x [0, Ly] of plate and marks which nodes belong to the plate.
6+% For plate.shape = 'rectangle' every node does. For 'ellipse' the plate is
7+% the ellipse inscribed in the box, (x - Lx/2)^2 / (Lx/2)^2 +
8+% (y - Ly/2)^2 / (Ly/2)^2 <= 1, a disk when Lx = Ly, and its nodes are those
9+% within half a grid spacing of the inside, so that the staircase boundary
10+% they form straddles the true edge rather than lying to one side of it.
11+% Nodes that do not belong to a grid cell with all four corners inside are
12+% dropped, so every node carries some area. Both finite-difference methods
13+% assemble their operators on this footprint (plate_matrices); the modal
14+% method has closed-form modes only for the rectangle.
15+%
16+% Fields.
17+% shape, nx, ny, dx, dy, x, y, X, Y the grid; w(i, j) lives at (x(j), y(i))
18+% inside ny x nx logical, the plate's nodes
19+% cell (ny-1) x (nx-1) logical, cells with all four corners inside
20+% m ny x nx trapezoidal node weights: adjacent plate cells / 4
21+% (1 inside, 1/2 on a straight edge, 1/4 at a convex corner)
22+% boundary inside nodes with a 4-neighbour outside: the ring where
23+% w = 0 is imposed for supported and clamped edges
24+% interior inside and not boundary
25+% full3 inside nodes whose eight neighbours are inside, where the
26+% curvatures can be taken by central differences
27+% near ny x nx index of the nearest full3 node, for every inside
28+% node (0 outside); plate_curvature_domain copies from it
29+% nunknowns_dirichlet, nunknowns_free sizes of the two unknown sets
30+% kappa, normal_x, normal_y ellipse only: curvature and outward unit
31+% normal of the true edge nearest each boundary node (zero
32+% elsewhere); the simply supported ellipse needs them, see
33+% plate_matrices
34+%
35+% Accuracy. A staircase approximates a curved edge to first order in the
36+% grid spacing: the geometry, the edge conditions and the stress near the
37+% edge all carry O(dx) errors there, against O(dx^2) inside. Refining the
38+% grid is the check.
39+
40+if ~isfield(plate, 'shape')
41+ shape = 'rectangle';
42+else
43+ shape = plate.shape;
44+end
45+if ~any(strcmp(shape, {'rectangle', 'ellipse'}))
46+ error('plate_domain: unknown plate shape "%s" (expected rectangle or ellipse)', shape);
47+end
48+Lx = plate.Lx; Ly = plate.Ly;
49+dx = Lx / (nx - 1);
50+dy = Ly / (ny - 1);
51+x = linspace(0, Lx, nx);
52+y = linspace(0, Ly, ny);
53+[X, Y] = meshgrid(x, y);
54+N = nx * ny;
55+
56+kappa = zeros(ny, nx);
57+normal_x = zeros(ny, nx);
58+normal_y = zeros(ny, nx);
59+if strcmp(shape, 'rectangle')
60+ in0 = true(ny, nx);
61+else
62+ xc = Lx / 2; yc = Ly / 2; a = Lx / 2; b = Ly / 2;
63+ phi = ((X - xc) / a).^2 + ((Y - yc) / b).^2; % 1 on the edge
64+ gx = 2 * (X - xc) / a^2; gy = 2 * (Y - yc) / b^2; % gradient of phi
65+ g = sqrt(gx.^2 + gy.^2);
66+ g(g == 0) = eps; % the centre
67+ d = (phi - 1) ./ g; % signed distance to the edge, first order
68+ hn = 1 ./ sqrt((gx ./ g / dx).^2 + (gy ./ g / dy).^2); % grid spacing along the normal
69+ in0 = d <= 0.5 * hn; % round the edge to the nearest node
70+end
71+cell = in0(1:end - 1, 1:end - 1) & in0(2:end, 1:end - 1) & in0(1:end - 1, 2:end) & in0(2:end, 2:end);
72+cnt = zeros(ny, nx);
73+cnt(1:end - 1, 1:end - 1) = cnt(1:end - 1, 1:end - 1) + cell;
74+cnt(2:end, 1:end - 1) = cnt(2:end, 1:end - 1) + cell;
75+cnt(1:end - 1, 2:end) = cnt(1:end - 1, 2:end) + cell;
76+cnt(2:end, 2:end) = cnt(2:end, 2:end) + cell;
77+inside = cnt > 0;
78+m = cnt / 4;
79+
80+% Neighbourhoods, with everything off the grid counting as outside
81+pad = false(ny + 2, nx + 2);
82+pad(2:end - 1, 2:end - 1) = inside;
83+nb4 = pad(1:end - 2, 2:end - 1) & pad(3:end, 2:end - 1) & pad(2:end - 1, 1:end - 2) & pad(2:end - 1, 3:end);
84+nb8 = nb4 & pad(1:end - 2, 1:end - 2) & pad(1:end - 2, 3:end) & pad(3:end, 1:end - 2) & pad(3:end, 3:end);
85+interior = inside & nb4;
86+boundary = inside & ~nb4;
87+full3 = inside & nb8;
88+if ~any(full3(:))
89+ error('plate_domain: the %d x %d grid is too coarse for this %s; no node has a full stencil inside the plate', nx, ny, shape);
90+end
91+
92+% Nearest node with a full stencil, by flooding outwards from those nodes
93+near = zeros(ny, nx);
94+kf = find(full3);
95+near(kf) = kf;
96+todo = inside & ~full3;
97+for pass = 1:max(nx, ny)
98+ if ~any(todo(:))
99+ break
100+ end
101+ np = zeros(ny + 2, nx + 2);
102+ np(2:end - 1, 2:end - 1) = near;
103+ shifts = [-1 0; 1 0; 0 -1; 0 1];
104+ for k = 1:4
105+ sh = np(2 + shifts(k, 1):end - 1 + shifts(k, 1), 2 + shifts(k, 2):end - 1 + shifts(k, 2));
106+ take = todo & sh > 0;
107+ near(take) = sh(take);
108+ todo(take) = false;
109+ end
110+end
111+
112+% Ellipse: curvature and normal of the true edge at the point nearest each
113+% boundary node (the radial projection of the node onto the ellipse)
114+if strcmp(shape, 'ellipse')
115+ kb = find(boundary);
116+ t = 1 ./ sqrt(max(phi(kb), eps));
117+ xe = xc + (X(kb) - xc) .* t;
118+ ye = yc + (Y(kb) - yc) .* t;
119+ px = 2 * (xe - xc) / a^2; py = 2 * (ye - yc) / b^2; % gradient of phi on the edge
120+ pg = sqrt(px.^2 + py.^2);
121+ normal_x(kb) = px ./ pg;
122+ normal_y(kb) = py ./ pg;
123+ kappa(kb) = ((2 / a^2) * py.^2 + (2 / b^2) * px.^2) ./ pg.^3; % curvature of the level set
124+end
125+
126+dom.shape = shape;
127+dom.nx = nx; dom.ny = ny; dom.dx = dx; dom.dy = dy;
128+dom.x = x; dom.y = y; dom.X = X; dom.Y = Y;
129+dom.inside = inside; dom.cell = cell; dom.m = m;
130+dom.boundary = boundary; dom.interior = interior; dom.full3 = full3; dom.near = near;
131+dom.nunknowns_dirichlet = nnz(interior);
132+dom.nunknowns_free = nnz(inside);
133+dom.kappa = kappa; dom.normal_x = normal_x; dom.normal_y = normal_y;
134+dom.N = N;
135+end
packages/core/src/m/plate_fd.mmodified+19−147View file
@@ -13,16 +13,27 @@ function result = plate_fd(params)
1313 % where w-, w, w+ are consecutive time levels and L is the five-point
1414 % Laplacian. Applying L twice gives the 13-point biharmonic stencil. The
1515 % scheme is second order in space and time (Bilbao, Numerical Sound
16-% Synthesis, 2009, ch. 12). Simply supported and clamped edges enter through
17-% one layer of ghost nodes (plate_pad). Free edges have no ghost-layer form,
18-% so for them L and L(L .) are the sparse energy-method operators of
19-% plate_matrices applied to every node, and the three rigid-body motions are
20-% projected out of the displacement after each step (plate_rigid_modes).
16+% Synthesis, 2009, ch. 12). On a rectangle, simply supported and clamped
17+% edges enter through one layer of ghost nodes (plate_pad). Free edges have
18+% no ghost-layer form, so for them L and L(L .) are the sparse energy-method
19+% operators of plate_matrices applied to every node, and the rigid-body
20+% motions are projected out of the displacement after each step
21+% (plate_rigid_modes). The same path serves the centre-pinned plate, whose
22+% pinned node is zeroed after each step, and the ellipse (plate_domain),
23+% whose edge is a staircase on the grid with every condition assembled as a
24+% sparse operator.
2125 %
2226 % The step is dt = cfl * dt_max with dt_max from plate_fd_plan; the scheme
2327 % is unstable for cfl > 1 and the run stops as soon as the field diverges,
2428 % reporting the step at which it did (info.diverged_step).
2529 %
30+% The work is in three stages, each its own file, so that a host can run the
31+% time loop elsewhere (the lab has a WebGPU port of plate_fd_loop) while the
32+% setup and the output stay this code:
33+% S = plate_fd_setup(params) grid, step, load, probe, listener, operators
34+% R = plate_fd_loop(S) the time stepping; records probe, frames, audio bins
35+% result = plate_fd_output(S, R) stresses, audio, spectra, the result struct
36+%
2637 % Approximations. Second-order truncation error in dx, dy, dt; the probe is
2738 % snapped to the nearest grid node; stress fields are evaluated only at the
2839 % saved frames, so the envelope is a sample of the true peak; the audio is
@@ -30,146 +41,7 @@ function result = plate_fd(params)
3041 %
3142 % See plate_result_schema for the fields of result.
3243
33-t_start = tic;
34-P = plate_setup(params);
35-plan = plate_fd_plan(params);
36-nx = plan.nx; ny = plan.ny; dx = plan.dx; dy = plan.dy; dt = plan.dt; nsteps = plan.nsteps;
37-x = linspace(0, P.plate.Lx, nx);
38-y = linspace(0, P.plate.Ly, ny);
39-[X, Y] = meshgrid(x, y); % w(i, j) lives at (x(j), y(i))
40-fs = P.output.fs;
41-notes = {};
42-if P.fd.cfl > 1
43- notes{end + 1} = sprintf('cfl = %g exceeds the stability limit; expect the solution to blow up', P.fd.cfl);
44-end
45-
46-s = plate_force_footprint(P.excitation, X, Y, dx * dy);
47-tt = (0:nsteps) * dt;
48-F = plate_force_time(P.excitation, tt);
49-
50-% Probe: nearest interior node
51-[~, jp] = min(abs(x - P.probe.x)); jp = min(max(jp, 2), nx - 1);
52-[~, ip] = min(abs(y - P.probe.y)); ip = min(max(ip, 2), ny - 1);
53-if abs(x(jp) - P.probe.x) > 0.5 * dx || abs(y(ip) - P.probe.y) > 0.5 * dy
54- notes{end + 1} = sprintf('probe moved to the nearest interior node (%.4g, %.4g) m', x(jp), y(ip));
55-end
56-
57-% Frames at the steps nearest to nframes equally spaced times in [0, duration]
58-frame_steps = unique(max(1, min(nsteps, round(linspace(0, P.fd.duration, P.output.nframes) / dt))));
59-nf = numel(frame_steps);
60-Wf = zeros(ny, nx, nf); VMf = zeros(ny, nx, nf); SPf = zeros(ny, nx, nf);
61-tf = zeros(1, nf);
62-
63-% Listener: every grid node is a patch; the acceleration of a step is added into
64-% the audio sample that contains it, delayed by the node's travel time.
65-L = plate_listener(P.listener, X, Y, dx * dy, fs);
66-bins = L.delay(:) + 1;
67-wgt = L.weight(:);
68-na = ceil(nsteps * dt * fs) + L.nbins + 1;
69-Pbins = zeros(L.nbins, na);
70-
71-% Free edges: sparse operators on every node, and the rigid-body projector
72-free = strcmp(P.bc, 'free');
73-if free
74- [Lm, Km, Mm] = plate_matrices(nx, ny, dx, dy, P.bc, P.material.nu);
75- mvec = full(diag(Mm));
76- [Rb, Qb] = plate_rigid_modes(X, Y, mvec);
77- notes{end + 1} = 'free edges: the rigid-body translation and tilts are removed from the displacement after every step';
78-end
79-
80-% State
81-w = zeros(ny, nx); wm = zeros(ny, nx); Lw_m = zeros(ny, nx);
82-w_p = zeros(1, nsteps + 1); wxx_p = zeros(1, nsteps + 1); wyy_p = zeros(1, nsteps + 1); wxy_p = zeros(1, nsteps + 1);
83-c_old = 1 - P.damping.sigma0 * dt;
84-c_new = 1 / (1 + P.damping.sigma0 * dt);
85-k2dt2 = P.kappa^2 * dt^2;
86-s1dt2 = 2 * P.damping.sigma1 * dt;
87-fdt2 = dt^2 / P.mu;
88-blowup = 1e3 * P.plate.h; % deflection beyond which the run is declared unstable
89-diverged_step = 0;
90-kf = 1;
91-
92-for nstep = 1:nsteps
93- if free
94- Lw = reshape((Lm * w(:)) ./ mvec, ny, nx);
95- B = reshape((Km * w(:)) ./ mvec, ny, nx);
96- else
97- Lw = plate_laplacian(plate_pad(w, P.bc), dx, dy);
98- B = zeros(ny, nx);
99- B(2:end - 1, 2:end - 1) = (Lw(3:end, 2:end - 1) - 2 * Lw(2:end - 1, 2:end - 1) + Lw(1:end - 2, 2:end - 1)) / dy^2 ...
100- + (Lw(2:end - 1, 3:end) - 2 * Lw(2:end - 1, 2:end - 1) + Lw(2:end - 1, 1:end - 2)) / dx^2;
101- end
102- wn = c_new * (2 * w - c_old * wm - k2dt2 * B + s1dt2 * (Lw - Lw_m) + fdt2 * F(nstep) * s);
103- if free
104- wn = wn - reshape(Rb * (Qb' * wn(:)), ny, nx);
105- else
106- wn(1, :) = 0; wn(end, :) = 0; wn(:, 1) = 0; wn(:, end) = 0;
107- end
108-
109- % Probe displacement and curvatures from the 3 x 3 neighbourhood
110- w_p(nstep + 1) = wn(ip, jp);
111- wxx_p(nstep + 1) = (wn(ip, jp + 1) - 2 * wn(ip, jp) + wn(ip, jp - 1)) / dx^2;
112- wyy_p(nstep + 1) = (wn(ip + 1, jp) - 2 * wn(ip, jp) + wn(ip - 1, jp)) / dy^2;
113- wxy_p(nstep + 1) = (wn(ip + 1, jp + 1) - wn(ip + 1, jp - 1) - wn(ip - 1, jp + 1) + wn(ip - 1, jp - 1)) / (4 * dx * dy);
114-
115- % Radiated pressure: acceleration at this step, weighted and delayed per node
116- acc = (wn - 2 * w + wm) / dt^2;
117- ia = floor(nstep * dt * fs) + 1;
118- Pbins(:, ia) = Pbins(:, ia) + accumarray(bins, wgt .* acc(:), [L.nbins 1]) * (dt * fs);
119-
120- if kf <= nf && nstep == frame_steps(kf)
121- if free
122- [wxx, wyy, wxy] = plate_curvature_free(wn, dx, dy, P.material.nu);
123- else
124- [wxx, wyy, wxy] = plate_curvature(plate_pad(wn, P.bc), dx, dy);
125- end
126- S = plate_stress(P, wxx, wyy, wxy);
127- Wf(:, :, kf) = wn; VMf(:, :, kf) = S.vm; SPf(:, :, kf) = S.sp; tf(kf) = nstep * dt;
128- kf = kf + 1;
129- end
130-
131- peak = max(abs(wn(:)));
132- if ~(peak < blowup) % also catches NaN and Inf
133- diverged_step = nstep;
134- notes{end + 1} = sprintf('solution diverged at step %d (t = %.3g s); results are truncated there', nstep, nstep * dt);
135- break
136- end
137- wm = w; w = wn; Lw_m = Lw;
138-end
139-
140-% Truncate after a blow-up so the outputs stay finite
141-if diverged_step > 0
142- nkeep = diverged_step; % keep the last finite level
143- tt = tt(1:nkeep); w_p = w_p(1:nkeep); wxx_p = wxx_p(1:nkeep); wyy_p = wyy_p(1:nkeep); wxy_p = wxy_p(1:nkeep);
144- keepf = frame_steps < diverged_step;
145- Wf = Wf(:, :, keepf); VMf = VMf(:, :, keepf); SPf = SPf(:, :, keepf); tf = tf(keepf);
146-end
147-
148-% Audio: sum the delay bins, then trim to the simulated duration
149-p_audio = zeros(1, na);
150-for d = 1:L.nbins
151- p_audio(d:end) = p_audio(d:end) + Pbins(d, 1:na - d + 1);
152-end
153-n_audio = min(na, ceil(numel(tt) * dt * fs));
154-p_audio = p_audio(1:n_audio);
155-
156-Sp = plate_stress(P, wxx_p, wyy_p, wxy_p);
157-[fw, Aw] = plate_spectrum(tt, w_p);
158-[~, Avm] = plate_spectrum(tt, Sp.vm);
159-
160-result.method = 'fd';
161-result.grid = struct('x', x, 'y', y);
162-result.t = tt;
163-result.probe = struct('x', x(jp), 'y', y(ip), 'w', w_p, 'sxx', Sp.sxx, 'syy', Sp.syy, 'txy', Sp.txy, 'vm', Sp.vm, 'sp', Sp.sp);
164-result.frames = struct('t', tf, 'w', Wf, 'vm', VMf, 'sp', SPf);
165-if isempty(tf)
166- result.envelope = struct('vm', zeros(ny, nx), 'sp', zeros(ny, nx));
167-else
168- result.envelope = struct('vm', max(VMf, [], 3), 'sp', max(abs(SPf), [], 3));
169-end
170-result.spectrum = struct('f', fw, 'w', Aw, 'vm', Avm);
171-result.audio = struct('fs', fs, 'p', p_audio, 't0', L.t0);
172-result.info = struct('nx', nx, 'ny', ny, 'dx', dx, 'dy', dy, 'dt', dt, 'dt_max', plan.dt_max, 'cfl', P.fd.cfl, ...
173- 'nsteps', nsteps, 'ncells', nx * ny, 'cell_steps', plan.cell_steps, 'fs', fs, ...
174- 'stable', diverged_step == 0, 'diverged_step', diverged_step, 'wall_s', toc(t_start), 'notes', {notes});
44+S = plate_fd_setup(params);
45+R = plate_fd_loop(S);
46+result = plate_fd_output(S, R);
17547 end
packages/core/src/m/plate_fd_loop.madded+111−0View file
@@ -0,0 +1,111 @@
1+function R = plate_fd_loop(S)
2+% PLATE_FD_LOOP The explicit time stepping of plate_fd.
3+%
4+% R = plate_fd_loop(S), with S from plate_fd_setup, advances
5+%
6+% wn = c_new * (2 w - c_old wm - k2dt2 B + s1dt2 (Lw - Lw_m) + fdt2 F(n) s)
7+%
8+% from rest for S.nsteps steps, where Lw is the five-point Laplacian of the
9+% current field, B the biharmonic (the Laplacian applied twice, the 13-point
10+% stencil) and Lw_m the Laplacian of the previous field. On a rectangle with
11+% supported or clamped edges the operators are the dense stencils with one
12+% layer of ghost nodes (plate_pad) and the edge nodes are held at zero.
13+% Otherwise (S.sparse) they are the sparse operators S.Lm, S.Km of
14+% plate_matrices on every node, divided by the node weights S.mvec; the
15+% nodes that are not unknowns (off the plate, on a supported ring, the pin)
16+% are zeroed with S.keep, and the S.nr rigid-body motions of a free plate
17+% are projected out of the new field with S.Rb, S.Qb.
18+%
19+% What the loop records, and returns in R:
20+% w_p, wxx_p, wyy_p, wxy_p displacement and curvatures at the probe node,
21+% 1 x (nsteps + 1); entry n + 1 is step n (t = n dt)
22+% Wf the field at the frame steps, ny x nx x nf
23+% Pbins radiated pressure per listener delay bin and
24+% audio sample, nbins x na: the acceleration of
25+% every node at every step, weighted by
26+% S.wgt and box-averaged onto audio samples
27+% diverged_step 0, or the step at which max |w| reached
28+% S.blowup (the run stops there)
29+% notes remarks about the run
30+% backend 'numbl': where the loop ran
31+% loop_s wall time of the loop by tic/toc
32+%
33+% The lab has a WebGPU port of this function (packages/core/src/gpu/fdLoop.ts)
34+% that computes the same record from the same S, in single precision. It
35+% carries the state as (w, w - wm) rather than (w, wm), which is the same
36+% arithmetic rearranged so that the rounding of w does not feed back into the
37+% velocity; keep the two implementations in step.
38+
39+t_loop = tic;
40+nx = S.nx; ny = S.ny; dx = S.dx; dy = S.dy; dt = S.dt; nsteps = S.nsteps; fs = S.fs;
41+bc = S.P.bc; use_sparse = S.sparse; nr = S.nr;
42+s = S.s; F = S.F; ip = S.ip; jp = S.jp;
43+frame_steps = S.frame_steps; nf = S.nf;
44+bins = S.bins; wgt = S.wgt; nbins = S.Lis.nbins; na = S.na;
45+c_old = S.c_old; c_new = S.c_new; k2dt2 = S.k2dt2; s1dt2 = S.s1dt2; fdt2 = S.fdt2; blowup = S.blowup;
46+if use_sparse
47+ Lm = S.Lm; Km = S.Km; mvec = S.mvec; keep = S.keep; Rb = S.Rb; Qb = S.Qb;
48+end
49+
50+% State
51+w = zeros(ny, nx); wm = zeros(ny, nx); Lw_m = zeros(ny, nx);
52+w_p = zeros(1, nsteps + 1); wxx_p = zeros(1, nsteps + 1); wyy_p = zeros(1, nsteps + 1); wxy_p = zeros(1, nsteps + 1);
53+Wf = zeros(ny, nx, nf);
54+Pbins = zeros(nbins, na);
55+diverged_step = 0;
56+notes = {};
57+kf = 1;
58+
59+for nstep = 1:nsteps
60+ if use_sparse
61+ Lw = reshape((Lm * w(:)) ./ mvec, ny, nx);
62+ B = reshape((Km * w(:)) ./ mvec, ny, nx);
63+ else
64+ Lw = plate_laplacian(plate_pad(w, bc), dx, dy);
65+ B = zeros(ny, nx);
66+ B(2:end - 1, 2:end - 1) = (Lw(3:end, 2:end - 1) - 2 * Lw(2:end - 1, 2:end - 1) + Lw(1:end - 2, 2:end - 1)) / dy^2 ...
67+ + (Lw(2:end - 1, 3:end) - 2 * Lw(2:end - 1, 2:end - 1) + Lw(2:end - 1, 1:end - 2)) / dx^2;
68+ end
69+ wn = c_new * (2 * w - c_old * wm - k2dt2 * B + s1dt2 * (Lw - Lw_m) + fdt2 * F(nstep) * s);
70+ if use_sparse
71+ wn = wn .* keep;
72+ if nr > 0
73+ wn = wn - reshape(Rb * (Qb' * wn(:)), ny, nx);
74+ end
75+ else
76+ wn(1, :) = 0; wn(end, :) = 0; wn(:, 1) = 0; wn(:, end) = 0;
77+ end
78+
79+ % Probe displacement and curvatures from the 3 x 3 neighbourhood
80+ w_p(nstep + 1) = wn(ip, jp);
81+ wxx_p(nstep + 1) = (wn(ip, jp + 1) - 2 * wn(ip, jp) + wn(ip, jp - 1)) / dx^2;
82+ wyy_p(nstep + 1) = (wn(ip + 1, jp) - 2 * wn(ip, jp) + wn(ip - 1, jp)) / dy^2;
83+ wxy_p(nstep + 1) = (wn(ip + 1, jp + 1) - wn(ip + 1, jp - 1) - wn(ip - 1, jp + 1) + wn(ip - 1, jp - 1)) / (4 * dx * dy);
84+
85+ % Radiated pressure: acceleration at this step, weighted and delayed per node
86+ acc = (wn - 2 * w + wm) / dt^2;
87+ ia = floor(nstep * dt * fs) + 1;
88+ Pbins(:, ia) = Pbins(:, ia) + accumarray(bins, wgt .* acc(:), [nbins 1]) * (dt * fs);
89+
90+ if kf <= nf && nstep == frame_steps(kf)
91+ Wf(:, :, kf) = wn;
92+ kf = kf + 1;
93+ end
94+
95+ peak = max(abs(wn(:)));
96+ if ~(peak < blowup) % also catches NaN and Inf
97+ diverged_step = nstep;
98+ notes{end + 1} = sprintf('solution diverged at step %d (t = %.3g s); results are truncated there', nstep, nstep * dt);
99+ break
100+ end
101+ wm = w; w = wn; Lw_m = Lw;
102+end
103+
104+R.w_p = w_p; R.wxx_p = wxx_p; R.wyy_p = wyy_p; R.wxy_p = wxy_p;
105+R.Wf = Wf;
106+R.Pbins = Pbins;
107+R.diverged_step = diverged_step;
108+R.notes = notes;
109+R.backend = 'numbl';
110+R.loop_s = toc(t_loop);
111+end
packages/core/src/m/plate_fd_output.madded+92−0View file
@@ -0,0 +1,92 @@
1+function result = plate_fd_output(S, R)
2+% PLATE_FD_OUTPUT The result struct of plate_fd from the record of its time loop.
3+%
4+% result = plate_fd_output(S, R), with S from plate_fd_setup and R from
5+% plate_fd_loop (or its WebGPU port), truncates the record after a blow-up,
6+% evaluates the surface stresses at the saved frames and at the probe, sums
7+% the listener's delay bins into the radiated pressure, takes the probe
8+% spectra and assembles the fields listed in plate_result_schema. The
9+% stresses are computed here, in double precision, from the displacement
10+% frames the loop kept, whichever precision the loop ran in. On an ellipse
11+% the nodes outside the plate are set to NaN in every field, so they are
12+% drawn blank.
13+
14+P = S.P;
15+dom = S.dom;
16+nx = S.nx; ny = S.ny; dx = S.dx; dy = S.dy; dt = S.dt; fs = S.fs;
17+tt = S.tt;
18+frame_steps = S.frame_steps;
19+w_p = R.w_p; wxx_p = R.wxx_p; wyy_p = R.wyy_p; wxy_p = R.wxy_p;
20+Wf = R.Wf;
21+notes = [S.notes, R.notes];
22+rect = strcmp(dom.shape, 'rectangle');
23+
24+% Truncate after a blow-up so the outputs stay finite
25+if R.diverged_step > 0
26+ nkeep = R.diverged_step; % keep the last finite level
27+ tt = tt(1:nkeep); w_p = w_p(1:nkeep); wxx_p = wxx_p(1:nkeep); wyy_p = wyy_p(1:nkeep); wxy_p = wxy_p(1:nkeep);
28+ keepf = frame_steps < R.diverged_step;
29+ Wf = Wf(:, :, keepf); frame_steps = frame_steps(keepf);
30+end
31+nf = numel(frame_steps);
32+tf = frame_steps * dt;
33+
34+% Stress at the saved frames, from the curvatures of each frame
35+VMf = zeros(ny, nx, nf); SPf = zeros(ny, nx, nf);
36+for kf = 1:nf
37+ if ~rect
38+ [wxx, wyy, wxy] = plate_curvature_domain(Wf(:, :, kf), dom);
39+ elseif S.free
40+ [wxx, wyy, wxy] = plate_curvature_free(Wf(:, :, kf), dx, dy, P.material.nu);
41+ else
42+ [wxx, wyy, wxy] = plate_curvature(plate_pad(Wf(:, :, kf), P.bc), dx, dy);
43+ end
44+ Sf = plate_stress(P, wxx, wyy, wxy);
45+ VMf(:, :, kf) = Sf.vm; SPf(:, :, kf) = Sf.sp;
46+end
47+if isempty(tf)
48+ env_vm = zeros(ny, nx); env_sp = zeros(ny, nx);
49+else
50+ env_vm = max(VMf, [], 3); env_sp = max(abs(SPf), [], 3);
51+end
52+if ~rect
53+ % Off the plate there is nothing to show
54+ out = ~dom.inside;
55+ for kf = 1:nf
56+ tmp = Wf(:, :, kf); tmp(out) = NaN; Wf(:, :, kf) = tmp;
57+ tmp = VMf(:, :, kf); tmp(out) = NaN; VMf(:, :, kf) = tmp;
58+ tmp = SPf(:, :, kf); tmp(out) = NaN; SPf(:, :, kf) = tmp;
59+ end
60+ env_vm(out) = NaN; env_sp(out) = NaN;
61+end
62+
63+% Audio: sum the delay bins, then trim to the simulated duration
64+na = S.na;
65+p_audio = zeros(1, na);
66+for d = 1:S.Lis.nbins
67+ p_audio(d:end) = p_audio(d:end) + R.Pbins(d, 1:na - d + 1);
68+end
69+n_audio = min(na, ceil(numel(tt) * dt * fs));
70+p_audio = p_audio(1:n_audio);
71+
72+Sp = plate_stress(P, wxx_p, wyy_p, wxy_p);
73+[fw, Aw] = plate_spectrum(tt, w_p);
74+[~, Avm] = plate_spectrum(tt, Sp.vm);
75+
76+result.method = 'fd';
77+result.grid = struct('x', S.x, 'y', S.y);
78+result.t = tt;
79+result.probe = struct('x', S.x(S.jp), 'y', S.y(S.ip), 'w', w_p, 'sxx', Sp.sxx, 'syy', Sp.syy, 'txy', Sp.txy, 'vm', Sp.vm, 'sp', Sp.sp);
80+result.frames = struct('t', tf, 'w', Wf, 'vm', VMf, 'sp', SPf);
81+result.envelope = struct('vm', env_vm, 'sp', env_sp);
82+result.spectrum = struct('f', fw, 'w', Aw, 'vm', Avm);
83+result.audio = struct('fs', fs, 'p', p_audio, 't0', S.Lis.t0);
84+result.info = struct('nx', nx, 'ny', ny, 'dx', dx, 'dy', dy, 'dt', dt, 'dt_max', S.plan.dt_max, 'cfl', P.fd.cfl, ...
85+ 'nsteps', S.nsteps, 'ncells', nx * ny, 'nunknowns', S.nunknowns, 'shape', dom.shape, 'cell_steps', S.plan.cell_steps, 'fs', fs, ...
86+ 'stable', R.diverged_step == 0, 'diverged_step', R.diverged_step, ...
87+ 'loop', R.backend, 'loop_s', R.loop_s, 'wall_s', toc(S.t_start), 'notes', {notes});
88+if ~isempty(S.fixed)
89+ result.info.support_x = S.X(S.fixed);
90+ result.info.support_y = S.Y(S.fixed);
91+end
92+end
packages/core/src/m/plate_fd_plan.mmodified+16−6View file
@@ -2,10 +2,11 @@ function plan = plate_fd_plan(params)
22 % PLATE_FD_PLAN Grid, time step and cost of the explicit finite-difference scheme.
33 %
44 % plan = plate_fd_plan(params) returns
5-% plan.nx, plan.ny, plan.dx, plan.dy grid (nodes include the edges)
5+% plan.nx, plan.ny, plan.dx, plan.dy grid (nodes include the edges of the bounding box)
66 % plan.dt_max largest stable time step [s]
77 % plan.dt step actually used, params.fd.cfl * dt_max
8-% plan.nsteps, plan.ncells, plan.cell_steps problem size
8+% plan.nsteps, plan.ncells, plan.cell_steps problem size (ncells counts every grid node)
9+% plan.nunknowns nodes the scheme actually updates
910 %
1011 % Stability. Von Neumann analysis of the scheme in plate_fd for a Fourier
1112 % mode on which the discrete Laplacian has eigenvalue -lambda gives the
@@ -19,10 +20,13 @@ function plan = plate_fd_plan(params)
1920 % Because dt scales with dx^2, halving the grid spacing costs eight times
2021 % as many cell updates for the same simulated time.
2122 %
22-% The same bound serves the free plate: the largest eigenvalue of its
23-% energy-method operator M \ K (plate_matrices) lies just below lambda_max^2
24-% (about 0.97 to 0.99 of it on the grids used here), so the step is safe and
25-% slightly conservative there.
23+% The same bound serves the free and centre-pinned plates and the ellipse:
24+% the Dirichlet Laplacian of any sub-domain has its spectrum inside that of
25+% the full grid, the clamped correction of plate_matrices does not raise the
26+% row sums above the interior stencil's, and the largest eigenvalue of the
27+% energy-method operator M \ K lies just below lambda_max^2 (about 0.97 to
28+% 0.99 of it on the grids used here), so the step is safe and slightly
29+% conservative there.
2630
2731 P = plate_setup(params);
2832 plan.nx = P.fd.nx;
@@ -38,5 +42,11 @@ plan.dt_max = 2 / (lambda_max * (s1 + sqrt(P.kappa^2 + s1^2)));
3842 plan.dt = P.fd.cfl * plan.dt_max;
3943 plan.nsteps = ceil(P.fd.duration / plan.dt);
4044 plan.ncells = plan.nx * plan.ny;
45+dom = plate_domain(P.plate, plan.nx, plan.ny);
46+if any(strcmp(P.bc, {'free', 'pinned-center'}))
47+ plan.nunknowns = dom.nunknowns_free;
48+else
49+ plan.nunknowns = dom.nunknowns_dirichlet;
50+end
4151 plan.cell_steps = plan.ncells * plan.nsteps;
4252 end
packages/core/src/m/plate_fd_setup.madded+152−0View file
@@ -0,0 +1,152 @@
1+function S = plate_fd_setup(params)
2+% PLATE_FD_SETUP Everything the explicit finite-difference time loop needs, computed once.
3+%
4+% S = plate_fd_setup(params) does the part of plate_fd that comes before the
5+% time stepping: the derived constants (plate_setup), the grid and time step
6+% (plate_fd_plan), the plate's footprint on the grid (plate_domain), the load
7+% footprint and force history, the probe node, the frame schedule, the
8+% listener geometry and, whenever the loop needs them, the sparse operators
9+% and the rigid-body projector (plate_matrices, plate_rigid_modes).
10+% plate_fd_loop advances the plate with S, and plate_fd_output turns what the
11+% loop recorded into the result struct; plate_fd calls the three in turn.
12+%
13+% The split lets a host run the loop somewhere else while the setup and the
14+% output stay this code: the lab has a WebGPU port of plate_fd_loop that reads
15+% exactly the fields below and hands its record to plate_fd_output.
16+%
17+% Two forms of the spatial operators. A rectangle with supported or clamped
18+% edges uses the dense stencils on the whole ny x nx array (plate_pad,
19+% plate_laplacian). Everything else, free edges on any shape and every
20+% condition on an ellipse, uses the sparse operators of plate_matrices on the
21+% unknowns, embedded in vectors over all N = nx ny nodes so that the loop
22+% keeps one array layout: rows and columns of nodes that are not unknowns are
23+% zero, their node weight is one, and S.keep zeroes them after every step.
24+%
25+% Fields.
26+% dom the footprint (plate_domain)
27+% nx, ny, dx, dy, x, y, X, Y grid; w(i, j) lives at (x(j), y(i))
28+% dt, nsteps, tt, F, fs time step, step count, tt = (0:nsteps) dt,
29+% force F(tt), audio sample rate
30+% s load footprint, ny x nx [1/m^2], zero off the plate
31+% ip, jp row and column of the probe node
32+% frame_steps, nf steps at which the field is kept
33+% Lis, bins, wgt, na plate_listener geometry; delay bin and
34+% weight of every node in column order; number
35+% of audio samples to accumulate
36+% free the edges are free (with or without the pin)
37+% nunknowns nodes in the system (the pinned node included)
38+% sparse the loop uses Lm, Km, mvec, keep, Rb, Qb
39+% Lm, Km, mvec N x N Laplacian and stiffness, N node weights
40+% keep ny x nx, 1 at the unknowns and 0 elsewhere
41+% (off the plate, on a supported ring, at the pin)
42+% nr, Rb, Qb rigid-body motions to remove (0, 2 or 3):
43+% N x nr basis and projector, zero off the plate
44+% fixed node index of the pinned node, or empty
45+% c_old, c_new, k2dt2, s1dt2, constants of the update
46+% fdt2, blowup wn = c_new (2 w - c_old wm - k2dt2 B
47+% + s1dt2 (Lw - Lw_m) + fdt2 F(n) s)
48+% with c_old = 1 - sigma0 dt, c_new =
49+% 1 / (1 + sigma0 dt), k2dt2 = kappa^2 dt^2,
50+% s1dt2 = 2 sigma1 dt, fdt2 = dt^2 / mu; blowup
51+% is the deflection that declares the run
52+% unstable
53+% notes remarks about this particular run
54+% t_start tic at entry, for the wall time
55+
56+S.t_start = tic;
57+P = plate_setup(params);
58+plan = plate_fd_plan(params);
59+nx = plan.nx; ny = plan.ny; dx = plan.dx; dy = plan.dy; dt = plan.dt; nsteps = plan.nsteps;
60+dom = plate_domain(P.plate, nx, ny);
61+x = dom.x; y = dom.y; X = dom.X; Y = dom.Y; % w(i, j) lives at (x(j), y(i))
62+N = nx * ny;
63+fs = P.output.fs;
64+notes = {};
65+if P.fd.cfl > 1
66+ notes{end + 1} = sprintf('cfl = %g exceeds the stability limit; expect the solution to blow up', P.fd.cfl);
67+end
68+rect = strcmp(dom.shape, 'rectangle');
69+free = any(strcmp(P.bc, {'free', 'pinned-center'}));
70+use_sparse = free || ~rect;
71+if free
72+ unknown = dom.inside;
73+else
74+ unknown = dom.interior;
75+end
76+
77+% Load on the unknowns only: a supported ring and the outside carry none
78+s = plate_force_footprint(P.excitation, X, Y, dx * dy, unknown);
79+tt = (0:nsteps) * dt;
80+F = plate_force_time(P.excitation, tt);
81+
82+% Probe: the nearest node whose eight neighbours belong to the plate, so the
83+% curvatures there are central differences
84+d2 = (X - P.probe.x).^2 + (Y - P.probe.y).^2;
85+d2(~dom.full3) = Inf;
86+[~, kp] = min(d2(:));
87+ip = mod(kp - 1, ny) + 1;
88+jp = floor((kp - 1) / ny) + 1;
89+if abs(x(jp) - P.probe.x) > 0.5 * dx || abs(y(ip) - P.probe.y) > 0.5 * dy
90+ notes{end + 1} = sprintf('probe moved to the nearest interior node (%.4g, %.4g) m', x(jp), y(ip));
91+end
92+
93+% Frames at the steps nearest to nframes equally spaced times in [0, duration]
94+frame_steps = unique(max(1, min(nsteps, round(linspace(0, P.fd.duration, P.output.nframes) / dt))));
95+
96+% Listener: every grid node is a patch; the acceleration of a step is added into
97+% the audio sample that contains it, delayed by the node's travel time.
98+Lis = plate_listener(P.listener, X, Y, dx * dy, fs);
99+
100+% Sparse operators on the unknowns, embedded over all N nodes, and the
101+% projector onto the complement of the rigid-body motions
102+fixed = zeros(0, 1);
103+Lm = []; Km = []; mvec = []; Rb = []; Qb = []; nr = 0;
104+keep = double(unknown);
105+nunknowns = nnz(unknown); % the pinned node counts: it is in the system, held at zero
106+if use_sparse
107+ [Lu, Ku, Mu, idx, fixed_u] = plate_matrices(dom, P.bc, P.material.nu);
108+ [ki, kj, kv] = find(Ku);
109+ Km = sparse(idx(ki), idx(kj), kv, N, N);
110+ [li, lj, lv] = find(Lu);
111+ Lm = sparse(idx(li), idx(lj), lv, N, N);
112+ mvec = ones(N, 1);
113+ mvec(idx) = full(diag(Mu));
114+ if ~isempty(fixed_u)
115+ fixed = idx(fixed_u);
116+ keep(fixed) = 0;
117+ end
118+ if free
119+ if isempty(fixed)
120+ [Ru, Qu] = plate_rigid_modes(X(idx), Y(idx), mvec(idx));
121+ notes{end + 1} = 'free edges: the rigid-body translation and tilts are removed from the displacement after every step';
122+ else
123+ [Ru, Qu] = plate_rigid_modes(X(idx), Y(idx), mvec(idx), [X(fixed) Y(fixed)]);
124+ notes{end + 1} = sprintf('pinned at the node (%.4g, %.4g) m nearest the centre; the two rigid rotations about it are removed after every step', X(fixed), Y(fixed));
125+ end
126+ nr = size(Ru, 2);
127+ Rb = zeros(N, nr); Rb(idx, :) = Ru;
128+ Qb = zeros(N, nr); Qb(idx, :) = Qu;
129+ end
130+end
131+if ~rect
132+ notes{end + 1} = sprintf('%s on a staircase of %d grid nodes; the edge is first order in the grid spacing', dom.shape, nnz(dom.inside));
133+end
134+
135+S.P = P; S.plan = plan; S.dom = dom;
136+S.nx = nx; S.ny = ny; S.dx = dx; S.dy = dy; S.x = x; S.y = y; S.X = X; S.Y = Y;
137+S.dt = dt; S.nsteps = nsteps; S.tt = tt; S.F = F; S.fs = fs;
138+S.s = s;
139+S.ip = ip; S.jp = jp;
140+S.frame_steps = frame_steps; S.nf = numel(frame_steps);
141+S.Lis = Lis; S.bins = Lis.delay(:) + 1; S.wgt = Lis.weight(:);
142+S.na = ceil(nsteps * dt * fs) + Lis.nbins + 1;
143+S.free = free; S.sparse = use_sparse; S.fixed = fixed; S.nunknowns = nunknowns;
144+S.Lm = Lm; S.Km = Km; S.mvec = mvec; S.keep = keep; S.nr = nr; S.Rb = Rb; S.Qb = Qb;
145+S.c_old = 1 - P.damping.sigma0 * dt;
146+S.c_new = 1 / (1 + P.damping.sigma0 * dt);
147+S.k2dt2 = P.kappa^2 * dt^2;
148+S.s1dt2 = 2 * P.damping.sigma1 * dt;
149+S.fdt2 = dt^2 / P.mu;
150+S.blowup = 1e3 * P.plate.h; % deflection beyond which the run is declared unstable
151+S.notes = notes;
152+end
packages/core/src/m/plate_force_footprint.mmodified+10−2View file
@@ -1,4 +1,4 @@
1-function s = plate_force_footprint(exc, X, Y, dA)
1+function s = plate_force_footprint(exc, X, Y, dA, mask)
22 % PLATE_FORCE_FOOTPRINT Spatial distribution s(x,y) [1/m^2] of the applied force.
33 %
44 % s = plate_force_footprint(exc, X, Y, dA) evaluates a raised-cosine bump of
@@ -6,6 +6,11 @@ function s = plate_force_footprint(exc, X, Y, dA)
66 % points X, Y (each with area dA) and scales it so that sum(s(:)) * dA = 1.
77 % The pressure is then f = F(t) s(x,y) and F(t) is the total force.
88 %
9+% s = plate_force_footprint(exc, X, Y, dA, mask) keeps only the points where
10+% the logical mask is true (the nodes of the plate that carry the load) and
11+% normalizes over those, so a strike near a curved edge or a supported edge
12+% still delivers the whole force.
13+%
914 % A finite footprint stands in for the contact patch of a hammer. It also
1015 % keeps the stress finite: a true point load gives a logarithmic stress
1116 % singularity in thin-plate theory, so the stress at the strike point would
@@ -14,10 +19,13 @@ function s = plate_force_footprint(exc, X, Y, dA)
1419 r = sqrt((X - exc.x).^2 + (Y - exc.y).^2);
1520 s = zeros(size(X));
1621 in = r < exc.radius;
22+if nargin >= 5
23+ in = in & mask;
24+end
1725 s(in) = 0.5 * (1 + cos(pi * r(in) / exc.radius));
1826 total = sum(s(:)) * dA;
1927 if total <= 0
20- error('plate_force_footprint: no grid point lies inside the footprint (radius %g m); enlarge the radius or refine the grid', exc.radius);
28+ error('plate_force_footprint: no grid point of the plate lies inside the footprint (radius %g m); enlarge the radius, refine the grid or move the strike inside the plate', exc.radius);
2129 end
2230 s = s / total;
2331 end
packages/core/src/m/plate_implicit.mmodified+63−34View file
@@ -11,19 +11,21 @@ function result = plate_implicit(params)
1111 % (w+ - 2w + w-)/dt^2 + sigma0 (w+ - w-)/dt - sigma1 L (w+ - w-)/dt
1212 % + kappa^2 B (theta w+ + (1 - 2 theta) w + theta w-) = F s / mu,
1313 %
14-% with L and B the sparse Laplacian and biharmonic matrices of plate_matrices,
15-% which also supplies the diagonal mass matrix M: the identity on the interior
16-% nodes for simply supported or clamped edges, the trapezoidal node weights on
17-% every node for free edges. Multiplying by dt^2 and collecting terms gives
18-% one linear system per step with a constant matrix,
14+% with L and B the sparse Laplacian and biharmonic matrices of plate_matrices
15+% on the plate's footprint (plate_domain), which also supplies the diagonal
16+% mass matrix M: the identity on the interior nodes for simply supported or
17+% clamped edges, the trapezoidal node weights on every node for free edges.
18+% Multiplying by dt^2 and collecting terms gives one linear system per step
19+% with a constant matrix,
1920 %
2021 % A w+ = C w - A_old w- + dt^2 F M s / mu
2122 % A = (1 + sigma0 dt) M - sigma1 dt L + theta kappa^2 dt^2 B
2223 % C = 2 M - (1 - 2 theta) kappa^2 dt^2 B
2324 % A_old = (1 - sigma0 dt) M + sigma1 dt L + theta kappa^2 dt^2 B.
2425 %
25-% For free edges the three rigid-body motions are projected out of w+ after
26-% every solve (plate_rigid_modes).
26+% For free edges the rigid-body motions are projected out of w+ after every
27+% solve (plate_rigid_modes); for the centre-pinned plate the pinned node is
28+% zeroed first.
2729 %
2830 % A is symmetric positive definite, so it is factored once by sparse Cholesky
2931 % and each step costs two sparse triangular solves plus two sparse products.
@@ -35,32 +37,30 @@ function result = plate_implicit(params)
3537 % Approximations. Second order in space and time. Implicit averaging
3638 % lengthens the period of a mode by about (omega dt)^2 / 12, so the highest
3739 % modes are detuned rather than unstable. Probe, frames, stress envelope and
38-% audio follow the conventions of plate_fd (see plate_result_schema).
40+% audio follow the conventions of plate_fd (see plate_result_schema); on an
41+% ellipse the edge is a staircase and first order, and the nodes outside the
42+% plate are NaN in the fields.
3943
4044 t_start = tic;
4145 P = plate_setup(params);
4246 plan = plate_implicit_plan(params);
4347 nx = plan.nx; ny = plan.ny; dx = plan.dx; dy = plan.dy; dt = plan.dt; nsteps = plan.nsteps;
4448 theta = P.implicit.theta;
45-x = linspace(0, P.plate.Lx, nx);
46-y = linspace(0, P.plate.Ly, ny);
47-[X, Y] = meshgrid(x, y);
49+dom = plate_domain(P.plate, nx, ny);
50+x = dom.x; y = dom.y; X = dom.X; Y = dom.Y;
4851 fs = P.output.fs;
4952 notes = {};
5053 if theta < 0.25
5154 notes{end + 1} = sprintf('theta = %g is below 1/4, so the scheme is only conditionally stable', theta);
5255 end
56+rect = strcmp(dom.shape, 'rectangle');
5357
5458 % ---- Matrices on the unknowns, factored once
55-free = strcmp(P.bc, 'free');
56-[L, B, M, idx, mx, my] = plate_matrices(nx, ny, dx, dy, P.bc, P.material.nu);
59+free = any(strcmp(P.bc, {'free', 'pinned-center'}));
60+[L, B, M, idx, fixed] = plate_matrices(dom, P.bc, P.material.nu);
5761 m = numel(idx);
5862 mvec = full(diag(M));
59-if free
60- rr = 1:ny; cc = 1:nx; % every node is an unknown
61-else
62- rr = 2:ny - 1; cc = 2:nx - 1; % the edge nodes carry w = 0
63-end
63+unknown = false(ny, nx); unknown(idx) = true;
6464 k2dt2 = P.kappa^2 * dt^2;
6565 A = (1 + P.damping.sigma0 * dt) * M - P.damping.sigma1 * dt * L + theta * k2dt2 * B;
6666 C = 2 * M - (1 - 2 * theta) * k2dt2 * B;
@@ -71,18 +71,27 @@ Rt = R'; % transposed once, not once per step
7171 factor_s = toc(t_factor);
7272
7373 % ---- Load, probe, frames, listener (as in plate_fd)
74-s = plate_force_footprint(P.excitation, X, Y, dx * dy);
75-s_unk = s(rr, cc);
76-load = (dt^2 / P.mu) * (mvec .* s_unk(:)); % M s on the unknowns
77-if free
78- [Rb, Qb] = plate_rigid_modes(X(rr, cc), Y(rr, cc), mvec);
74+s = plate_force_footprint(P.excitation, X, Y, dx * dy, unknown);
75+load = (dt^2 / P.mu) * (mvec .* s(idx)); % M s on the unknowns
76+if free && isempty(fixed)
77+ [Rb, Qb] = plate_rigid_modes(X(idx), Y(idx), mvec);
7978 notes{end + 1} = 'free edges: the rigid-body translation and tilts are removed from the displacement after every step';
79+elseif free
80+ [Rb, Qb] = plate_rigid_modes(X(idx), Y(idx), mvec, [X(idx(fixed)) Y(idx(fixed))]);
81+ notes{end + 1} = sprintf('pinned at the node (%.4g, %.4g) m nearest the centre; the two rigid rotations about it are removed after every step', X(idx(fixed)), Y(idx(fixed)));
82+end
83+if ~rect
84+ notes{end + 1} = sprintf('%s on a staircase of %d grid nodes; the edge is first order in the grid spacing', dom.shape, nnz(dom.inside));
8085 end
8186 tt = (0:nsteps) * dt;
8287 F = plate_force_time(P.excitation, tt);
8388
84-[~, jp] = min(abs(x - P.probe.x)); jp = min(max(jp, 2), nx - 1);
85-[~, ip] = min(abs(y - P.probe.y)); ip = min(max(ip, 2), ny - 1);
89+% Probe: the nearest node whose eight neighbours belong to the plate
90+d2 = (X - P.probe.x).^2 + (Y - P.probe.y).^2;
91+d2(~dom.full3) = Inf;
92+[~, kp] = min(d2(:));
93+ip = mod(kp - 1, ny) + 1;
94+jp = floor((kp - 1) / ny) + 1;
8695 if abs(x(jp) - P.probe.x) > 0.5 * dx || abs(y(ip) - P.probe.y) > 0.5 * dy
8796 notes{end + 1} = sprintf('probe moved to the nearest interior node (%.4g, %.4g) m', x(jp), y(ip));
8897 end
@@ -110,22 +119,27 @@ for nstep = 1:nsteps
110119 rhs = C * w - A_old * wm + F(nstep) * load;
111120 wn = R \ (Rt \ rhs); % forward then back substitution
112121 if free
122+ if ~isempty(fixed)
123+ wn(fixed) = 0;
124+ end
113125 wn = wn - Rb * (Qb' * wn);
114126 end
115127
116- wfull(rr, cc) = reshape(wn, my, mx);
128+ wfull(idx) = wn;
117129 w_p(nstep + 1) = wfull(ip, jp);
118130 wxx_p(nstep + 1) = (wfull(ip, jp + 1) - 2 * wfull(ip, jp) + wfull(ip, jp - 1)) / dx^2;
119131 wyy_p(nstep + 1) = (wfull(ip + 1, jp) - 2 * wfull(ip, jp) + wfull(ip - 1, jp)) / dy^2;
120132 wxy_p(nstep + 1) = (wfull(ip + 1, jp + 1) - wfull(ip + 1, jp - 1) - wfull(ip - 1, jp + 1) + wfull(ip - 1, jp - 1)) / (4 * dx * dy);
121133
122134 acc = zeros(ny, nx);
123- acc(rr, cc) = reshape((wn - 2 * w + wm) / dt^2, my, mx);
135+ acc(idx) = (wn - 2 * w + wm) / dt^2;
124136 ia = floor(nstep * dt * fs) + 1;
125137 Pbins(:, ia) = Pbins(:, ia) + accumarray(bins, wgt .* acc(:), [Lis.nbins 1]) * (dt * fs);
126138
127139 if kf <= nf && nstep == frame_steps(kf)
128- if free
140+ if ~rect
141+ [wxx, wyy, wxy] = plate_curvature_domain(wfull, dom);
142+ elseif free
129143 [wxx, wyy, wxy] = plate_curvature_free(wfull, dx, dy, P.material.nu);
130144 else
131145 [wxx, wyy, wxy] = plate_curvature(plate_pad(wfull, P.bc), dx, dy);
@@ -150,6 +164,21 @@ if diverged_step > 0
150164 keepf = frame_steps < diverged_step;
151165 Wf = Wf(:, :, keepf); VMf = VMf(:, :, keepf); SPf = SPf(:, :, keepf); tf = tf(keepf);
152166 end
167+nf = numel(tf);
168+if isempty(tf)
169+ env_vm = zeros(ny, nx); env_sp = zeros(ny, nx);
170+else
171+ env_vm = max(VMf, [], 3); env_sp = max(abs(SPf), [], 3);
172+end
173+if ~rect
174+ out = ~dom.inside;
175+ for kf = 1:nf
176+ tmp = Wf(:, :, kf); tmp(out) = NaN; Wf(:, :, kf) = tmp;
177+ tmp = VMf(:, :, kf); tmp(out) = NaN; VMf(:, :, kf) = tmp;
178+ tmp = SPf(:, :, kf); tmp(out) = NaN; SPf(:, :, kf) = tmp;
179+ end
180+ env_vm(out) = NaN; env_sp(out) = NaN;
181+end
153182
154183 p_audio = zeros(1, na);
155184 for d = 1:Lis.nbins
@@ -167,15 +196,15 @@ result.grid = struct('x', x, 'y', y);
167196 result.t = tt;
168197 result.probe = struct('x', x(jp), 'y', y(ip), 'w', w_p, 'sxx', Sp.sxx, 'syy', Sp.syy, 'txy', Sp.txy, 'vm', Sp.vm, 'sp', Sp.sp);
169198 result.frames = struct('t', tf, 'w', Wf, 'vm', VMf, 'sp', SPf);
170-if isempty(tf)
171- result.envelope = struct('vm', zeros(ny, nx), 'sp', zeros(ny, nx));
172-else
173- result.envelope = struct('vm', max(VMf, [], 3), 'sp', max(abs(SPf), [], 3));
174-end
199+result.envelope = struct('vm', env_vm, 'sp', env_sp);
175200 result.spectrum = struct('f', fw, 'w', Aw, 'vm', Avm);
176201 result.audio = struct('fs', fs, 'p', p_audio, 't0', Lis.t0);
177202 result.info = struct('nx', nx, 'ny', ny, 'dx', dx, 'dy', dy, 'dt', dt, 'theta', theta, ...
178- 'steps_per_sample', P.implicit.steps_per_sample, 'nunknowns', m, 'nnz_factor', nnz(R), 'factor_s', factor_s, ...
203+ 'steps_per_sample', P.implicit.steps_per_sample, 'nunknowns', m, 'shape', dom.shape, 'nnz_factor', nnz(R), 'factor_s', factor_s, ...
179204 'nsteps', nsteps, 'ncells', nx * ny, 'cell_steps', plan.cell_steps, 'fs', fs, ...
180205 'stable', diverged_step == 0, 'diverged_step', diverged_step, 'wall_s', toc(t_start), 'notes', {notes});
206+if ~isempty(fixed)
207+ result.info.support_x = X(idx(fixed));
208+ result.info.support_y = Y(idx(fixed));
209+end
181210 end
packages/core/src/m/plate_implicit_plan.mmodified+4−3View file
@@ -21,10 +21,11 @@ plan.dy = P.plate.Ly / (plan.ny - 1);
2121 plan.dt = 1 / (P.output.fs * P.implicit.steps_per_sample);
2222 plan.nsteps = ceil(P.implicit.duration / plan.dt);
2323 plan.ncells = plan.nx * plan.ny;
24-if strcmp(P.bc, 'free')
25- plan.nunknowns = plan.nx * plan.ny; % every node moves; see plate_matrices
24+dom = plate_domain(P.plate, plan.nx, plan.ny);
25+if any(strcmp(P.bc, {'free', 'pinned-center'}))
26+ plan.nunknowns = dom.nunknowns_free; % every node of the plate is in the system; see plate_matrices
2627 else
27- plan.nunknowns = (plan.nx - 2) * (plan.ny - 2);
28+ plan.nunknowns = dom.nunknowns_dirichlet; % the boundary ring carries w = 0
2829 end
2930 plan.cell_steps = plan.ncells * plan.nsteps;
3031 plan.theta = P.implicit.theta;
packages/core/src/m/plate_matrices.mmodified+202−83View file
@@ -1,102 +1,221 @@
1-function [L, K, M, idx, mx, my] = plate_matrices(nx, ny, dx, dy, bc, nu)
1+function [L, K, M, idx, fixed] = plate_matrices(dom, bc, nu)
22 % PLATE_MATRICES Sparse Laplacian, stiffness and mass matrices for the plate.
33 %
4-% [L, K, M, idx, mx, my] = plate_matrices(nx, ny, dx, dy, bc, nu) assembles
5-% the spatial operators of the finite-difference schemes as sparse matrices
6-% acting on the vector w(idx) of unknowns, where idx indexes the ny x nx field
7-% with y running fastest (column-major, as w(:) does) and the unknowns form an
8-% my x mx array. The semi-discrete equation is
4+% [L, K, M, idx, fixed] = plate_matrices(dom, bc, nu) assembles the spatial
5+% operators of the finite-difference schemes on the footprint dom from
6+% plate_domain, as sparse matrices acting on the vector w(idx) of unknowns.
7+% idx lists the unknown nodes in column-major order of the ny x nx field
8+% (y running fastest, as w(:) does). fixed is the unknown held at w = 0
9+% anyway (empty except for the pinned centre); the solvers zero it after
10+% every step. The semi-discrete equation is
911 %
1012 % M w_tt + 2 sigma0 M w_t - 2 sigma1 L w_t + kappa^2 K w = F(t) M s / mu
1113 %
12-% Simply supported and clamped edges (w = 0 on the edge nodes). The unknowns
13-% are the mx = nx-2 by my = ny-2 interior nodes, M is the identity, L the
14-% five-point Laplacian and K the thirteen-point biharmonic, with the edges
15-% handled exactly as in the ghost-node scheme of plate_fd / plate_pad: D2*D2
16-% encodes the odd reflection (w_nn = 0), and for clamped edges the even
17-% reflection changes the end coefficient of the fourth difference from 5/dx^4
18-% to 7/dx^4. The mixed term 2 w_xxyy never reaches a ghost node. nu is not
19-% used. A test in the lab checks K against the explicit stencil.
14+% Simply supported and clamped edges (w = 0 on the boundary ring). The
15+% unknowns are dom.interior, M is the identity, L the five-point Laplacian
16+% with zero Dirichlet data and K = L L, the thirteen-point biharmonic. On
17+% the rectangle this is exactly the ghost-node scheme of plate_fd and
18+% plate_pad: L L encodes the odd reflection (w_nn = 0). For clamped edges
19+% the even reflection makes the Laplacian nonzero on the ring, and the term
2020 %
21-% Free edges. Every node is an unknown (mx = nx, my = ny). The free-edge
22-% conditions, no bending moment and no Kirchhoff shear on an edge and no
23-% twisting moment at a corner, involve third derivatives and do not fit in a
24-% ghost layer. K is therefore assembled from the bending energy
21+% K = L L + sum_q m_q v_q v_q'
22+%
23+% adds it: v_q is the Laplacian at ring node q with each outside neighbour
24+% mirrored onto the inside one, and m_q the trapezoidal weight of q. On the
25+% rectangle this is the 7/dx^4 end coefficient in place of 5/dx^4; on a
26+% staircase edge it treats each face of the staircase as clamped.
27+%
28+% Simply supported curved edges and the plate paradox. With w = 0 on the
29+% boundary the bending energy is
30+%
31+% U = D/2 * [ integral (lap w)^2 dA - (1 - nu) * contour kappa w_n^2 ds ]
32+%
33+% where kappa is the curvature of the edge. L L is the first term alone.
34+% On a polygon, and on a staircase, kappa vanishes, so the discrete plate
35+% converges to one with lap w = 0 on the edge instead of zero bending
36+% moment: the Babuska (Sapondzhyan) paradox, which puts the fundamental of a
37+% simply supported disk 16 percent too high. For the ellipse the missing
38+% term is added with the curvature and normal of the true edge nearest each
39+% ring node (dom.kappa), w_n taken by one-sided differences into the plate
40+% and the arc length shared out over the faces of the staircase. It is a
41+% negative correction, small against the entries of L L but of order one
42+% in the energy of a smooth mode, and it restores the simply supported
43+% limit. Nothing is added for the rectangle, whose edges are straight.
44+%
45+% Free edges. Every node of the plate is an unknown (dom.inside). The
46+% free-edge conditions, no bending moment and no Kirchhoff shear on an edge
47+% and no twisting moment at a corner, involve third derivatives and do not
48+% fit in a ghost layer. K is therefore assembled from the bending energy
2549 %
2650 % U = D/2 * integral( w_xx^2 + w_yy^2 + 2 nu w_xx w_yy + 2 (1-nu) w_xy^2 ) dA
2751 %
28-% using second differences at the nodes, the cross difference at the cell
29-% centres and trapezoidal quadrature weights (the finite-difference energy
30-% method; Bushnell, Computers & Structures 1973). The free-edge conditions
31-% are then natural boundary conditions of the discrete energy: the scheme is
32-% second order inside and first order at the boundary. M holds the trapezoidal
33-% node weights (1 inside, 1/2 on an edge, 1/4 at a corner) and L is the
34-% Laplacian energy operator with the natural zero-normal-derivative condition.
35-% K is symmetric positive semidefinite and its null space is exactly the three
36-% rigid-body motions 1, x and y, which plate_rigid_modes removes from the
37-% state. All three matrices are divided by dx*dy, so inside the plate M = 1
38-% and K and L reduce to the usual stencils, and the largest eigenvalue of
39-% M \ K stays below the interior bound (4/dx^2 + 4/dy^2)^2 that plate_fd_plan
52+% using second differences at the nodes that have both neighbours, the
53+% cross difference at the cell centres and the trapezoidal node weights
54+% dom.m (the finite-difference energy method; Bushnell, Computers &
55+% Structures 1973). The free-edge conditions are then natural boundary
56+% conditions of the discrete energy: the scheme is second order inside and
57+% first order at the boundary. M holds the node weights and L is the
58+% Laplacian energy operator with the natural zero-normal-derivative
59+% condition, each grid edge weighted by its share of cells. K is symmetric
60+% positive semidefinite and its null space is exactly the three rigid-body
61+% motions 1, x and y, which plate_rigid_modes removes from the state. All
62+% three matrices are per unit cell area, so inside the plate M = 1 and K
63+% and L reduce to the usual stencils, and the largest eigenvalue of M \ K
64+% stays near the interior bound (4/dx^2 + 4/dy^2)^2 that plate_fd_plan
4065 % uses for the time step.
4166 %
42-% Everything is built from Kronecker products so the matrices stay sparse;
43-% kron(Ax, By) acts on w(:) as By * W * Ax'.
67+% Pinned centre (a cymbal on a stand). The edges are free and the unknown
68+% nearest the centre is held at w = 0: its row and column are zeroed in K
69+% and L, so its equation decouples (M keeps its unit diagonal and A stays
70+% positive definite) and the other nodes never see it. One pinned point
71+% leaves the two rigid rotations about it at zero frequency;
72+% plate_rigid_modes with the pin as pivot removes them, standing in for the
73+% felt washers of a stand. Under a point support the curvature grows like
74+% log r, so the stress right at the support depends on the grid, as it does
75+% under a point load.
4476
45-if nargin < 6
77+if nargin < 3
4678 nu = NaN;
4779 end
48-switch bc
49- case {'simply-supported', 'clamped'}
50- mx = nx - 2;
51- my = ny - 2;
52- ex = ones(mx, 1);
53- ey = ones(my, 1);
54- D2x = spdiags([ex -2 * ex ex], -1:1, mx, mx) / dx^2;
55- D2y = spdiags([ey -2 * ey ey], -1:1, my, my) / dy^2;
56- Ix = speye(mx);
57- Iy = speye(my);
58- D4x = D2x * D2x;
59- D4y = D2y * D2y;
60- if strcmp(bc, 'clamped')
61- D4x(1, 1) = 7 / dx^4; D4x(mx, mx) = 7 / dx^4;
62- D4y(1, 1) = 7 / dy^4; D4y(my, my) = 7 / dy^4;
80+nx = dom.nx; ny = dom.ny; dx = dom.dx; dy = dom.dy;
81+dirichlet = any(strcmp(bc, {'simply-supported', 'clamped'}));
82+energy = any(strcmp(bc, {'free', 'pinned-center'}));
83+if ~dirichlet && ~energy
84+ error('plate_matrices: unknown boundary condition "%s"', bc);
85+end
86+if energy && ~isfinite(nu)
87+ error('plate_matrices: free edges need Poisson''s ratio nu');
88+end
89+if dirichlet
90+ unk = dom.interior;
91+else
92+ unk = dom.inside;
93+end
94+idx = find(unk);
95+n = numel(idx);
96+me = (1:n)';
97+
98+% Unknown number of every node (0 for none), on a padded grid so that every
99+% node has four neighbours to look up
100+num_in = zeros(ny, nx);
101+num_in(idx) = me;
102+num = zeros(ny + 2, nx + 2);
103+num(2:end - 1, 2:end - 1) = num_in;
104+nbL = num(2:end - 1, 1:end - 2); nbR = num(2:end - 1, 3:end); % neighbours at x - dx, x + dx
105+nbD = num(1:end - 2, 2:end - 1); nbU = num(3:end, 2:end - 1); % neighbours at y - dy, y + dy
106+uL = nbL(idx); uR = nbR(idx); uD = nbD(idx); uU = nbU(idx); % of the unknowns, in unknown order
107+
108+if dirichlet
109+ % ---- Dirichlet Laplacian: full diagonal, off-diagonals only to unknowns
110+ ii = me; jj = me; vv = -(2 / dx^2 + 2 / dy^2) * ones(n, 1);
111+ nbrs = {uL, uR, uD, uU};
112+ coef = [1 / dx^2, 1 / dx^2, 1 / dy^2, 1 / dy^2];
113+ for k = 1:4
114+ nb = nbrs{k};
115+ has = find(nb > 0);
116+ ii = [ii; me(has)]; jj = [jj; nb(has)]; vv = [vv; coef(k) * ones(numel(has), 1)];
117+ end
118+ L = sparse(ii, jj, vv, n, n);
119+ K = L * L;
120+ M = speye(n);
121+ fixed = zeros(0, 1);
122+
123+ % Ring nodes and their unknown neighbours
124+ kb = find(dom.boundary);
125+ nb = numel(kb);
126+ qb = (1:nb)';
127+ bL = nbL(kb); bR = nbR(kb); bD = nbD(kb); bU = nbU(kb);
128+ if strcmp(bc, 'clamped')
129+ % Laplacian at each ring node with the outside neighbour mirrored inside
130+ ii = []; jj = []; vv = [];
131+ for k = 1:2
132+ if k == 1
133+ p1 = bL; p2 = bR; c = 1 / dx^2;
134+ else
135+ p1 = bD; p2 = bU; c = 1 / dy^2;
136+ end
137+ both = find(p1 > 0 & p2 > 0);
138+ one1 = find(p1 > 0 & p2 == 0);
139+ one2 = find(p2 > 0 & p1 == 0);
140+ ii = [ii; qb(both); qb(both); qb(one1); qb(one2)];
141+ jj = [jj; p1(both); p2(both); p1(one1); p2(one2)];
142+ vv = [vv; c * ones(2 * numel(both), 1); 2 * c * ones(numel(one1) + numel(one2), 1)];
63143 end
64- L = kron(Ix, D2y) + kron(D2x, Iy);
65- K = kron(Ix, D4y) + 2 * kron(D2x, D2y) + kron(D4x, Iy);
66- M = speye(mx * my);
67- mask = false(ny, nx);
68- mask(2:end - 1, 2:end - 1) = true;
69- idx = find(mask);
70- case 'free'
144+ V = sparse(ii, jj, vv, nb, n);
145+ mq = dom.m(kb);
146+ K = K + V' * spdiags(mq, 0, nb, nb) * V;
147+ elseif any(dom.kappa(:) ~= 0)
148+ % Simply supported curved edge: - (1 - nu) contour kappa w_n^2 ds, with
149+ % w_n^2 = w_x^2 + w_y^2 on an edge where w = 0, the gradient by one-sided
150+ % (or, with unknowns on both sides, central) differences at each ring
151+ % node, and ds the faces of the staircase projected onto the true edge
71152 if ~isfinite(nu)
72- error('plate_matrices: free edges need Poisson''s ratio nu');
153+ error('plate_matrices: a simply supported curved edge needs Poisson''s ratio nu');
154+ end
155+ ds = dy * abs(dom.normal_x(kb)) .* ((bL > 0) + (bR > 0)) + dx * abs(dom.normal_y(kb)) .* ((bD > 0) + (bU > 0));
156+ wq = (1 - nu) * dom.kappa(kb) .* ds / (dx * dy);
157+ for k = 1:2
158+ if k == 1
159+ p1 = bL; p2 = bR; h = dx;
160+ else
161+ p1 = bD; p2 = bU; h = dy;
162+ end
163+ both = find(p1 > 0 & p2 > 0);
164+ one1 = find(p1 > 0 & p2 == 0);
165+ one2 = find(p2 > 0 & p1 == 0);
166+ ii = [qb(both); qb(both); qb(one1); qb(one2)];
167+ jj = [p1(both); p2(both); p1(one1); p2(one2)];
168+ vv = [-ones(numel(both), 1) / (2 * h); ones(numel(both), 1) / (2 * h); ones(numel(one1) + numel(one2), 1) / h];
169+ G = sparse(ii, jj, vv, nb, n);
170+ K = K - G' * spdiags(wq, 0, nb, nb) * G;
73171 end
74- mx = nx;
75- my = ny;
76- % One-dimensional pieces, as triplets so the layout is unambiguous:
77- % second differences at the interior nodes ((n-2) x n), first differences
78- % at the cell midpoints ((n-1) x n), selectors of the interior nodes, and
79- % trapezoidal node weights.
80- kx = (1:nx - 2)'; ky = (1:ny - 2)';
81- d2x = sparse([kx; kx; kx], [kx; kx + 1; kx + 2], [ones(nx - 2, 1); -2 * ones(nx - 2, 1); ones(nx - 2, 1)], nx - 2, nx) / dx^2;
82- d2y = sparse([ky; ky; ky], [ky; ky + 1; ky + 2], [ones(ny - 2, 1); -2 * ones(ny - 2, 1); ones(ny - 2, 1)], ny - 2, ny) / dy^2;
83- jx = (1:nx - 1)'; jy = (1:ny - 1)';
84- d1x = sparse([jx; jx], [jx; jx + 1], [-ones(nx - 1, 1); ones(nx - 1, 1)], nx - 1, nx) / dx;
85- d1y = sparse([jy; jy], [jy; jy + 1], [-ones(ny - 1, 1); ones(ny - 1, 1)], ny - 1, ny) / dy;
86- Sx = sparse(kx, kx + 1, ones(nx - 2, 1), nx - 2, nx);
87- Sy = sparse(ky, ky + 1, ones(ny - 2, 1), ny - 2, ny);
88- tx = ones(nx, 1); tx([1 end]) = 0.5; Tx = spdiags(tx, 0, nx, nx);
89- ty = ones(ny, 1); ty([1 end]) = 0.5; Ty = spdiags(ty, 0, ny, ny);
90- % The four terms of the energy density, each as D' * (weights) * D
91- Kxx = kron(d2x' * d2x, Ty); % w_xx^2 at nodes 2..nx-1, every row
92- Kyy = kron(Tx, d2y' * d2y); % w_yy^2 at nodes 2..ny-1, every column
93- Kc = nu * (kron(d2x' * Sx, Sy' * d2y) + kron(Sx' * d2x, d2y' * Sy)); % 2 nu w_xx w_yy where both exist
94- Kxy = 2 * (1 - nu) * kron(d1x' * d1x, d1y' * d1y); % 2 (1-nu) w_xy^2 at the cell centres
95- K = Kxx + Kyy + Kc + Kxy;
96- L = -(kron(d1x' * d1x, Ty) + kron(Tx, d1y' * d1y));
97- M = kron(Tx, Ty);
98- idx = (1:nx * ny)';
99- otherwise
100- error('plate_matrices: unknown boundary condition "%s"', bc);
172+ end
173+else
174+ % ---- Energy method on every node of the plate
175+ mu = dom.m(idx);
176+ hasX = uL > 0 & uR > 0;
177+ hasY = uD > 0 & uU > 0;
178+ kx = find(hasX); nkx = numel(kx);
179+ ky = find(hasY); nky = numel(ky);
180+ Dxx = sparse([(1:nkx)'; (1:nkx)'; (1:nkx)'], [uL(kx); me(kx); uR(kx)], [ones(nkx, 1); -2 * ones(nkx, 1); ones(nkx, 1)] / dx^2, nkx, n);
181+ Dyy = sparse([(1:nky)'; (1:nky)'; (1:nky)'], [uD(ky); me(ky); uU(ky)], [ones(nky, 1); -2 * ones(nky, 1); ones(nky, 1)] / dy^2, nky, n);
182+ Kxx = Dxx' * spdiags(mu(kx), 0, nkx, nkx) * Dxx; % w_xx^2 where both x-neighbours exist
183+ Kyy = Dyy' * spdiags(mu(ky), 0, nky, nky) * Dyy; % w_yy^2 where both y-neighbours exist
184+ kc = find(hasX & hasY); nkc = numel(kc);
185+ Dxc = sparse([(1:nkc)'; (1:nkc)'; (1:nkc)'], [uL(kc); me(kc); uR(kc)], [ones(nkc, 1); -2 * ones(nkc, 1); ones(nkc, 1)] / dx^2, nkc, n);
186+ Dyc = sparse([(1:nkc)'; (1:nkc)'; (1:nkc)'], [uD(kc); me(kc); uU(kc)], [ones(nkc, 1); -2 * ones(nkc, 1); ones(nkc, 1)] / dy^2, nkc, n);
187+ Wc = spdiags(mu(kc), 0, nkc, nkc);
188+ Kc = nu * (Dxc' * Wc * Dyc + Dyc' * Wc * Dxc); % 2 nu w_xx w_yy where both exist
189+ % Cells with four unknown corners: the twist at the cell centre
190+ A00 = num_in(1:end - 1, 1:end - 1); A10 = num_in(2:end, 1:end - 1);
191+ A01 = num_in(1:end - 1, 2:end); A11 = num_in(2:end, 2:end);
192+ cu = A00 > 0 & A10 > 0 & A01 > 0 & A11 > 0;
193+ kcell = find(cu); ncell = numel(kcell);
194+ Dxy = sparse([(1:ncell)'; (1:ncell)'; (1:ncell)'; (1:ncell)'], [A00(kcell); A10(kcell); A01(kcell); A11(kcell)], ...
195+ [ones(ncell, 1); -ones(ncell, 1); -ones(ncell, 1); ones(ncell, 1)] / (dx * dy), ncell, n);
196+ Kxy = 2 * (1 - nu) * (Dxy' * Dxy); % 2 (1-nu) w_xy^2 at the cell centres
197+ K = Kxx + Kyy + Kc + Kxy;
198+ % Laplacian energy: first differences along the grid edges, each weighted
199+ % by its share of adjacent cells (1 inside, 1/2 along the boundary)
200+ cpad = zeros(ny + 1, nx + 1);
201+ cpad(2:end - 1, 2:end - 1) = cu;
202+ WX = (cpad(2:end, 2:end) + cpad(1:end - 1, 2:end)) / 2; % x-edge from (i,j) to (i,j+1)
203+ WY = (cpad(2:end, 2:end) + cpad(2:end, 1:end - 1)) / 2; % y-edge from (i,j) to (i+1,j)
204+ ex = find(uR > 0); nex = numel(ex);
205+ ey = find(uU > 0); ney = numel(ey);
206+ Dx1 = sparse([(1:nex)'; (1:nex)'], [me(ex); uR(ex)], [-ones(nex, 1); ones(nex, 1)] / dx, nex, n);
207+ Dy1 = sparse([(1:ney)'; (1:ney)'], [me(ey); uU(ey)], [-ones(ney, 1); ones(ney, 1)] / dy, ney, n);
208+ wx = WX(idx(ex)); wy = WY(idx(ey));
209+ L = -(Dx1' * spdiags(wx, 0, nex, nex) * Dx1 + Dy1' * spdiags(wy, 0, ney, ney) * Dy1);
210+ M = spdiags(mu, 0, n, n);
211+ fixed = zeros(0, 1);
212+ if strcmp(bc, 'pinned-center')
213+ xc = dom.x(end) / 2; yc = dom.y(end) / 2;
214+ [~, fixed] = min((dom.X(idx) - xc).^2 + (dom.Y(idx) - yc).^2); % the unknown nearest the centre
215+ sel = ones(n, 1); sel(fixed) = 0;
216+ Pz = spdiags(sel, 0, n, n);
217+ K = Pz * K * Pz; % its row and column drop out
218+ L = Pz * L * Pz;
219+ end
101220 end
102221 end
packages/core/src/m/plate_modal.mmodified+6−0View file
@@ -23,6 +23,9 @@ function result = plate_modal(params)
2323 % second derivatives of the modes, and the radiated pressure from the modal
2424 % accelerations through the Rayleigh integral (plate_listener).
2525 %
26+% The rectangle is the only shape with closed-form modes for these edges;
27+% an ellipse is refused, and the finite-difference methods handle it.
28+%
2629 % Approximations. The modal sum is truncated at params.modal.fmax and
2730 % params.modal.max_modes; the force is sampled at params.output.fs; the
2831 % footprint integral and the Rayleigh integral use quadrature grids. Stress
@@ -36,6 +39,9 @@ P = plate_setup(params);
3639 if ~strcmp(P.bc, 'simply-supported')
3740 error('plate_modal: closed-form modes exist only for simply supported edges (got "%s"); use the finite-difference method for other edges', P.bc);
3841 end
42+if ~strcmp(P.plate.shape, 'rectangle')
43+ error('plate_modal: closed-form modes exist only for the rectangle (got "%s"); use a finite-difference method for the ellipse', P.plate.shape);
44+end
3945 a = P.plate.Lx;
4046 b = P.plate.Ly;
4147 fs = P.output.fs;
packages/core/src/m/plate_result_schema.mmodified+7−4View file
@@ -1,18 +1,21 @@
11 % PLATE_RESULT_SCHEMA Fields shared by every solver's result struct.
22 %
3-% result.method 'modal' or 'fd'
4-% result.grid x (1 x nx), y (1 x ny): display grid of the frames [m]
3+% result.method 'modal', 'fd' or 'implicit'
4+% result.grid x (1 x nx), y (1 x ny): display grid of the frames [m],
5+% spanning the bounding box of the plate
56 % result.t probe time axis [s]
67 % result.probe x, y (actual probe location used), and the series
78 % w [m], sxx, syy, txy, vm, sp [Pa] sampled on result.t
89 % result.frames t (1 x nf) and the fields w [m], vm, sp [Pa],
9-% each ny x nx x nf (rows index y, columns index x)
10+% each ny x nx x nf (rows index y, columns index x);
11+% on an ellipse the nodes outside the plate are NaN
1012 % result.envelope vm, sp: peak over the saved frames, ny x nx [Pa]
1113 % result.spectrum f [Hz], w [m], vm [Pa]: amplitude spectra of the probe series
1214 % result.audio fs [Hz], p (pressure at the listener) [Pa], t0 (dropped travel time) [s]
1315 % result.info problem size, time step, wall_s (solve wall time, tic/toc inside numbl),
1416 % stable, diverged_step, notes (cell array of strings about approximations
15-% that were applied to this particular run)
17+% that were applied to this particular run); the finite-difference
18+% methods add shape and nunknowns
1619 % result.modes (modal only) m, n, f_hz, sigma, G, truncated
1720 %
1821 % This file documents the convention and has no executable content.
packages/core/src/m/plate_rigid_modes.mmodified+13−4View file
@@ -1,4 +1,4 @@
1-function [R, Q] = plate_rigid_modes(X, Y, mvec)
1+function [R, Q] = plate_rigid_modes(X, Y, mvec, pivot)
22 % PLATE_RIGID_MODES Basis of the rigid-body motions of a free plate and its M-orthogonal projector.
33 %
44 % [R, Q] = plate_rigid_modes(X, Y, mvec) takes the node coordinates X, Y and
@@ -12,6 +12,11 @@ function [R, Q] = plate_rigid_modes(X, Y, mvec)
1212 % product: what remains when the mean displacement and the two mean tilts are
1313 % taken out.
1414 %
15+% [R, Q] = plate_rigid_modes(X, Y, mvec, pivot), with pivot = [x0 y0] the
16+% coordinates of a pinned node, keeps only the two rotations about that point,
17+% R = [x - x0, y - y0]. Both columns vanish at the pin, so the projection
18+% leaves w = 0 there untouched.
19+%
1520 % Why remove them. A completely free plate has three zero-frequency modes,
1621 % and a strike gives it net momentum and angular momentum: with only the loss
1722 % sigma0 to stop it, the plate drifts to a displacement of order
@@ -26,9 +31,13 @@ function [R, Q] = plate_rigid_modes(X, Y, mvec)
2631 x = X(:);
2732 y = Y(:);
2833 m = mvec(:);
29-xc = sum(m .* x) / sum(m);
30-yc = sum(m .* y) / sum(m);
31-R = [ones(numel(x), 1), x - xc, y - yc];
34+if nargin >= 4
35+ R = [x - pivot(1), y - pivot(2)];
36+else
37+ xc = sum(m .* x) / sum(m);
38+ yc = sum(m .* y) / sum(m);
39+ R = [ones(numel(x), 1), x - xc, y - yc];
40+end
3241 G = R' * (m .* R); % 3 x 3 Gram matrix in the mass inner product
3342 Q = (m .* R) / G;
3443 end
packages/core/src/m/plate_setup.mmodified+13−2View file
@@ -8,6 +8,10 @@ function P = plate_setup(params)
88 % P.Kb bending-stress factor at the surface z = h/2 [N/m]
99 % P.Ks shear-stress factor at the surface z = h/2 [N/m]
1010 %
11+% The plate is params.plate.shape, 'rectangle' or 'ellipse', filling the
12+% box [0, Lx] x [0, Ly]: for the ellipse Lx and Ly are the diameters along
13+% x and y, and Lx = Ly is a disk. A missing shape means a rectangle.
14+%
1115 % The transverse displacement w(x,y,t) of the plate obeys
1216 %
1317 % w_tt + 2 sigma0 w_t - 2 sigma1 lap(w_t) + kappa^2 lap(lap(w)) = f(x,y,t) / mu
@@ -21,14 +25,21 @@ P = params;
2125 pl = params.plate;
2226 mat = params.material;
2327
28+if ~isfield(pl, 'shape')
29+ pl.shape = 'rectangle';
30+ P.plate.shape = 'rectangle';
31+end
32+if ~any(strcmp(pl.shape, {'rectangle', 'ellipse'}))
33+ error('plate_setup: unknown plate shape "%s" (expected rectangle or ellipse)', pl.shape);
34+end
2435 if pl.Lx <= 0 || pl.Ly <= 0 || pl.h <= 0
2536 error('plate_setup: plate dimensions must be positive');
2637 end
2738 if pl.h > 0.1 * min(pl.Lx, pl.Ly)
2839 warning('plate_setup: thickness is more than a tenth of the shorter side; Kirchhoff-Love thin-plate theory is not reliable here');
2940 end
30-if ~any(strcmp(params.bc, {'simply-supported', 'clamped', 'free'}))
31- error('plate_setup: unknown boundary condition "%s" (expected simply-supported, clamped or free)', params.bc);
41+if ~any(strcmp(params.bc, {'simply-supported', 'clamped', 'free', 'pinned-center'}))
42+ error('plate_setup: unknown boundary condition "%s" (expected simply-supported, clamped, free or pinned-center)', params.bc);
3243 end
3344
3445 P.D = mat.E * pl.h^3 / (12 * (1 - mat.nu^2));
packages/core/src/methods/common.tsmodified+23−0View file
@@ -33,6 +33,29 @@ export const FREE_EDGE_EQUATIONS: Equation[] = [
3333 },
3434 ];
3535
36+/** Free edges with the centre held, as a cymbal sits on its stand. */
37+export const PINNED_CENTER_EQUATIONS: Equation[] = [
38+ {
39+ tex: String.raw`w(x_c, y_c, t) = 0, \qquad \text{free edges}, \qquad \mathbf K \leftarrow \mathsf P \mathbf K \mathsf P`,
40+ caption:
41+ "Pinned centre: the node nearest the centre is held at zero displacement by zeroing its row and column of the assembled stiffness (P is the selector of the other nodes), so its equation decouples. One pinned point still leaves the two rigid rotations about it at zero frequency; they are removed after every step, which is the role the felt washers of a stand play. Under a point support the curvature grows like log r, so the stress right at the support depends on the grid, as under a point load.",
42+ },
43+];
44+
45+/** How the two finite-difference methods handle an elliptical plate on the rectangular grid. */
46+export const CURVED_EDGE_EQUATIONS: Equation[] = [
47+ {
48+ tex: String.raw`\Omega_h = \left\{ (x_j, y_i) : \frac{(x_j - L_x/2)^2}{(L_x/2)^2} + \frac{(y_i - L_y/2)^2}{(L_y/2)^2} \le 1 + O(\Delta x) \right\}`,
49+ caption:
50+ "Ellipse: the plate is the set of grid nodes within half a spacing of the inscribed ellipse, so its edge is a staircase that straddles the true curve. The same stencils and the same sparse assembly apply on this footprint; supported and clamped edges hold the outermost ring of nodes at w = 0, free edges use the energy method on every node. Geometry, edge conditions and the stress along the edge are then first order in the spacing (second order inside); refining the grid is the check.",
51+ },
52+ {
53+ tex: String.raw`U = \frac{D}{2}\iint (\nabla^2 w)^2\, dA \;-\; \frac{D(1-\nu)}{2}\oint \kappa\, w_n^2\, ds \qquad (w = 0 \text{ on the edge})`,
54+ caption:
55+ "Simply supported curved edge. With w = 0 on the edge the bending energy takes this form, kappa being the curvature of the edge. On a polygon or a staircase kappa vanishes, and a scheme that imposes only w = 0 converges to a plate with zero Laplacian on the edge instead of zero moment: the Babuska (Sapondzhyan) paradox, which puts the fundamental of a simply supported disk 16 percent too high. The lab adds the boundary term with the curvature of the true ellipse at each edge node, which restores the simply supported limit; for the rectangle the term is zero.",
56+ },
57+];
58+
3659 export const STRESS_EQUATIONS: Equation[] = [
3760 {
3861 tex: String.raw`\sigma_{xx} = -\frac{E z}{1-\nu^2}\left(w_{xx} + \nu w_{yy}\right),\quad \sigma_{yy} = -\frac{E z}{1-\nu^2}\left(w_{yy} + \nu w_{xx}\right),\quad \tau_{xy} = -\frac{E z}{1+\nu}\,w_{xy}`,
packages/core/src/methods/fd.tsmodified+23−4View file
@@ -1,7 +1,9 @@
11 import type { MethodDefinition } from "./types";
22 import {
3+ CURVED_EDGE_EQUATIONS,
34 FREE_EDGE_EQUATIONS,
45 GOVERNING_EQUATIONS,
6+ PINNED_CENTER_EQUATIONS,
57 RADIATION_EQUATIONS,
68 STRESS_EQUATIONS,
79 } from "./common";
@@ -14,20 +16,26 @@ export const FD_METHOD: MethodDefinition = {
1416 solverFunction: "plate_fd",
1517 planFunction: "plate_fd_plan",
1618 supports: {
17- bc: ["simply-supported", "clamped", "free"],
19+ shape: ["rectangle", "ellipse"],
20+ bc: ["simply-supported", "clamped", "free", "pinned-center"],
1821 excitation: ["impulse", "harmonic"],
1922 audio: "clip",
2023 },
2124 readingOrder: [
2225 "plate_fd.m",
26+ "plate_fd_setup.m",
27+ "plate_fd_loop.m",
28+ "plate_fd_output.m",
2329 "plate_fd_plan.m",
2430 "plate_setup.m",
31+ "plate_domain.m",
2532 "plate_pad.m",
2633 "plate_laplacian.m",
2734 "plate_curvature.m",
2835 "plate_matrices.m",
2936 "plate_rigid_modes.m",
3037 "plate_curvature_free.m",
38+ "plate_curvature_domain.m",
3139 "plate_force_time.m",
3240 "plate_force_footprint.m",
3341 "plate_stress.m",
@@ -45,27 +53,38 @@ export const FD_METHOD: MethodDefinition = {
4553 {
4654 tex: String.raw`\mathsf L w_{ij} = \frac{w_{i,j+1} - 2w_{ij} + w_{i,j-1}}{\Delta x^2} + \frac{w_{i+1,j} - 2w_{ij} + w_{i-1,j}}{\Delta y^2}`,
4755 caption:
48- "Supported and clamped edges enter through one ghost layer: odd reflection for simply supported (w = w_nn = 0), even reflection for clamped (w = w_n = 0). Free edges take the energy-method operators below.",
56+ "On a rectangle, supported and clamped edges enter through one ghost layer: odd reflection for simply supported (w = w_nn = 0), even reflection for clamped (w = w_n = 0). Free edges, and every condition on an ellipse, take the sparse operators below.",
4957 },
5058 {
5159 tex: String.raw`\Delta t \le \frac{2}{\lambda_{\max}\left(\sigma_1 + \sqrt{\kappa^2 + \sigma_1^2}\right)}, \qquad \lambda_{\max} = \frac{4}{\Delta x^2} + \frac{4}{\Delta y^2}`,
5260 caption:
53- "Stability bound from von Neumann analysis; for a square grid without sigma1 it reads dt <= dx^2 / (4 kappa). Because dt scales with dx^2, halving the spacing costs eight times more work.",
61+ "Stability bound from von Neumann analysis; for a square grid without sigma1 it reads dt <= dx^2 / (4 kappa). Because dt scales with dx^2, halving the spacing costs eight times more work. The sparse operators of the free plate and of the ellipse stay just under the same bound.",
5462 },
5563 ...FREE_EDGE_EQUATIONS,
64+ ...PINNED_CENTER_EQUATIONS,
65+ ...CURVED_EDGE_EQUATIONS,
5666 ...STRESS_EQUATIONS,
5767 ...RADIATION_EQUATIONS,
5868 ],
5969 validity: [
6070 "Second-order truncation error in dx, dy and dt; the top modes on the grid are dispersive and their frequencies are too low.",
61- "The probe is snapped to the nearest interior grid node.",
71+ "The probe is snapped to the nearest grid node with a full stencil inside the plate.",
6272 "Stress fields are evaluated only at the saved frames, so the envelope samples the true peak.",
6373 "The time-step limit makes long simulated times expensive: cost grows as (1/dx)^4 for a fixed duration.",
6474 "Audio is a short clip covering only the simulated window, box-averaged onto audio samples.",
75+ "With the GPU option the time loop runs in single precision (a WebGPU port of plate_fd_loop.m); the setup, the stresses, the audio summation and the spectra stay in double precision in numbl. Against the double-precision loop the displacement differs by a relative L2 error of about 1e-6 over 50 ms on a 61 x 41 grid, the stress and the sound by about 1e-4.",
6576 "Free edges: the moment and shear conditions are natural conditions of a discrete bending energy, met to first order at the boundary and second order inside; the three rigid-body motions are removed after every step, so the displacement is that of a plate on soft supports, not one drifting in space.",
77+ "Pinned centre: the support is one grid node, so the stress right at it is resolution-dependent, and the two rigid rotations about it are removed rather than restrained by a modelled washer.",
78+ "Ellipse: the edge is a staircase of grid nodes, so the geometry, the edge conditions and the stress along the edge are first order in the grid spacing. The stress on the outermost ring of nodes is copied from the nearest node with a full stencil, which puts the peak stress of a clamped edge about one spacing inside. For simply supported edges the curvature term that the staircase loses is added back (the plate paradox), so the disk's frequencies converge to the exact ones.",
6679 ],
6780 durationOf: (spec) => spec.fd.duration,
6881 gridOf: (spec) => [spec.fd.nx, spec.fd.ny],
6982 settingsKey: "fd",
7083 usPerCellStep: 0.35,
84+ gpu: {
85+ loopFile: "plate_fd_loop.m",
86+ setupFunction: "plate_fd_setup",
87+ outputFunction: "plate_fd_output",
88+ usPerCellStep: 0.02,
89+ },
7190 };
packages/core/src/methods/implicit.tsmodified+12−3View file
@@ -1,7 +1,9 @@
11 import type { MethodDefinition } from "./types";
22 import {
3+ CURVED_EDGE_EQUATIONS,
34 FREE_EDGE_EQUATIONS,
45 GOVERNING_EQUATIONS,
6+ PINNED_CENTER_EQUATIONS,
57 RADIATION_EQUATIONS,
68 STRESS_EQUATIONS,
79 } from "./common";
@@ -15,16 +17,19 @@ export const IMPLICIT_METHOD: MethodDefinition = {
1517 solverFunction: "plate_implicit",
1618 planFunction: "plate_implicit_plan",
1719 supports: {
18- bc: ["simply-supported", "clamped", "free"],
20+ shape: ["rectangle", "ellipse"],
21+ bc: ["simply-supported", "clamped", "free", "pinned-center"],
1922 excitation: ["impulse", "harmonic"],
2023 audio: "clip",
2124 },
2225 readingOrder: [
2326 "plate_implicit.m",
2427 "plate_implicit_plan.m",
28+ "plate_domain.m",
2529 "plate_matrices.m",
2630 "plate_rigid_modes.m",
2731 "plate_curvature_free.m",
32+ "plate_curvature_domain.m",
2833 "plate_setup.m",
2934 "plate_force_time.m",
3035 "plate_force_footprint.m",
@@ -40,7 +45,7 @@ export const IMPLICIT_METHOD: MethodDefinition = {
4045 {
4146 tex: String.raw`\frac{w^{n+1} - 2w^n + w^{n-1}}{\Delta t^2} + \sigma_0\frac{w^{n+1} - w^{n-1}}{\Delta t} - \sigma_1 \mathsf L\frac{w^{n+1} - w^{n-1}}{\Delta t} + \kappa^2 \mathsf B\left(\theta w^{n+1} + (1-2\theta) w^n + \theta w^{n-1}\right) = \frac{F^n s}{\mu}`,
4247 caption:
43- "Theta scheme: the biharmonic term is averaged over three time levels. L and B are the sparse Laplacian and biharmonic matrices built from the same stencils as the explicit scheme.",
48+ "Theta scheme: the biharmonic term is averaged over three time levels. L and B are the sparse Laplacian and biharmonic matrices built from the same stencils as the explicit scheme, on the plate's footprint.",
4449 },
4550 {
4651 tex: String.raw`\mathbf A\, w^{n+1} = \mathbf C\, w^n - \mathbf A_{\mathrm{old}}\, w^{n-1} + \frac{\Delta t^2}{\mu} F^n s, \qquad \mathbf A = (1+\sigma_0\Delta t)\mathbf I - \sigma_1\Delta t\,\mathsf L + \theta\kappa^2\Delta t^2\,\mathsf B`,
@@ -53,16 +58,20 @@ export const IMPLICIT_METHOD: MethodDefinition = {
5358 "No stability limit, so the step is chosen for accuracy. The price is period elongation: modes near half the sample rate run slow by tens of percent unless several steps per sample are used.",
5459 },
5560 ...FREE_EDGE_EQUATIONS,
61+ ...PINNED_CENTER_EQUATIONS,
62+ ...CURVED_EDGE_EQUATIONS,
5663 ...STRESS_EQUATIONS,
5764 ...RADIATION_EQUATIONS,
5865 ],
5966 validity: [
6067 "Second-order truncation error in dx and dy, identical to the explicit scheme, since the spatial stencils are the same.",
6168 "Second order in dt but dispersive: high modes are detuned rather than unstable. Two steps per sample keeps a 10 kHz mode within about 4 percent.",
62- "The probe is snapped to the nearest interior grid node; stress fields are evaluated only at the saved frames.",
69+ "The probe is snapped to the nearest grid node with a full stencil inside the plate; stress fields are evaluated only at the saved frames.",
6370 "Cost per step is a sparse triangular solve pair, about a millisecond for a 61 x 41 grid in the browser, independent of the plate stiffness.",
6471 "Audio covers only the simulated window.",
6572 "Free edges: the moment and shear conditions are natural conditions of a discrete bending energy, met to first order at the boundary and second order inside; the three rigid-body motions are removed after every step, so the displacement is that of a plate on soft supports, not one drifting in space.",
73+ "Pinned centre: the support is one grid node, so the stress right at it is resolution-dependent, and the two rigid rotations about it are removed rather than restrained by a modelled washer.",
74+ "Ellipse: the edge is a staircase of grid nodes, so the geometry, the edge conditions and the stress along the edge are first order in the grid spacing. The stress on the outermost ring of nodes is copied from the nearest node with a full stencil, which puts the peak stress of a clamped edge about one spacing inside. For simply supported edges the curvature term that the staircase loses is added back (the plate paradox), so the disk's frequencies converge to the exact ones.",
6675 ],
6776 durationOf: (spec) => spec.implicit.duration,
6877 gridOf: (spec) => [spec.implicit.nx, spec.implicit.ny],
packages/core/src/methods/modal.tsmodified+7−2View file
@@ -7,7 +7,12 @@ export const MODAL_METHOD: MethodDefinition = {
77 shortName: "Modal",
88 summary: "Closed-form modes of the simply supported rectangle, each integrated exactly in time.",
99 solverFunction: "plate_modal",
10- supports: { bc: ["simply-supported"], excitation: ["impulse", "harmonic"], audio: "full" },
10+ supports: {
11+ shape: ["rectangle"],
12+ bc: ["simply-supported"],
13+ excitation: ["impulse", "harmonic"],
14+ audio: "full",
15+ },
1116 readingOrder: [
1217 "plate_modal.m",
1318 "plate_setup.m",
@@ -40,7 +45,7 @@ export const MODAL_METHOD: MethodDefinition = {
4045 ...RADIATION_EQUATIONS,
4146 ],
4247 validity: [
43- "Exact only for simply supported edges; other edge conditions have no closed-form modes.",
48+ "Exact only for the simply supported rectangle; other edge conditions and the ellipse have no closed-form modes of this kind (a disk would need Bessel functions).",
4449 "The modal sum is truncated at fmax and max_modes. Displacement converges quickly with mode count; stress converges more slowly, most of all near the strike point.",
4550 "The footprint projection and the Rayleigh integral are evaluated by quadrature on fixed grids.",
4651 "Damping rates are a model with two parameters, chosen, not measured.",
packages/core/src/methods/registry.tsmodified+5−1View file
@@ -15,9 +15,13 @@ export function getMethod(id: string): MethodDefinition {
1515 /** Reasons a method cannot run on a spec (empty when it can). */
1616 export function methodBlockers(method: MethodDefinition, spec: ProblemSpec): string[] {
1717 const out: string[] = [];
18+ if (!method.supports.shape.includes(spec.plate.shape))
19+ out.push(
20+ `${method.name} has no solution for the ${spec.plate.shape}; the finite-difference methods embed it in their grid.`,
21+ );
1822 if (!method.supports.bc.includes(spec.bc))
1923 out.push(
20- `${method.name} does not support ${spec.bc} edges (supports: ${method.supports.bc.join(", ")}).`,
24+ `${method.name} does not support the edge condition "${spec.bc}" (supports: ${method.supports.bc.join(", ")}).`,
2125 );
2226 if (!method.supports.excitation.includes(spec.excitation.type))
2327 out.push(`${method.name} does not support the ${spec.excitation.type} excitation.`);
packages/core/src/methods/types.tsmodified+14−1View file
@@ -1,4 +1,4 @@
1-import type { BoundaryCondition, ExcitationType, ProblemSpec } from "../spec/types";
1+import type { BoundaryCondition, ExcitationType, PlateShape, ProblemSpec } from "../spec/types";
22
33 export interface Equation {
44 /** LaTeX, display mode */
@@ -25,6 +25,7 @@ export interface MethodDefinition {
2525 /** Optional `.m` function taking `params` and returning a cost plan */
2626 planFunction?: string;
2727 supports: {
28+ shape: PlateShape[];
2829 bc: BoundaryCondition[];
2930 excitation: ExcitationType[];
3031 audio: AudioSupport;
@@ -42,4 +43,16 @@ export interface MethodDefinition {
4243 settingsKey: "modal" | "fd" | "implicit";
4344 /** Starting guess for wall time per cell-step in microseconds, refined by measurements */
4445 usPerCellStep: number;
46+ /**
47+ * Set when the method's time loop has a WebGPU port: the loop file it mirrors
48+ * (an edit there sends the run back to numbl) and the setup and output
49+ * functions the host runs in numbl around it.
50+ */
51+ gpu?: {
52+ loopFile: string;
53+ setupFunction: string;
54+ outputFunction: string;
55+ /** Starting guess for the GPU path, per cell-step, refined by measurements */
56+ usPerCellStep: number;
57+ };
4558 }
packages/core/src/node.tsmodified+2−0View file
@@ -1,3 +1,5 @@
11 /** Node-only entry: the in-process runner used by tests and scripts. */
22 export { NodeRunner } from "./runner/nodeRunner";
33 export { executeBundle, numblVersion } from "./runner/execute";
4+export { executeRequest, type GpuProvider } from "./runner/executeGpu";
5+export { requestGpuDevice, runFdLoopGpu, type FdLoopSetup, type FdLoopRecord } from "./gpu/fdLoop";
packages/core/src/runner/executeGpu.tsadded+121−0View file
@@ -0,0 +1,121 @@
1+import type { BackendInfo, ComputeInfo, RunRequest, RunResult } from "./types";
2+import { executeBundle, formatError } from "./execute";
3+import { LIBRARY_FILES } from "../generate/library";
4+import { runFdLoopGpu, type GpuHandle } from "../gpu/fdLoop";
5+import { fdLoopRecordToPlain, fdLoopSetupFromPlain } from "../gpu/fdStage";
6+
7+/** Supplies the (cached) WebGPU device of the host, or rejects when there is none. */
8+export type GpuProvider = () => Promise<GpuHandle>;
9+
10+/**
11+ * Run a request, on the GPU when it asks for that and the host has one. The
12+ * plain path is executeBundle. The staged path runs the setup script in numbl,
13+ * the loop on the GPU and the output script in numbl, and reports where the
14+ * loop ran in `compute`; when the GPU was asked for but cannot be used, the
15+ * whole bundle runs in numbl and `compute.note` says why.
16+ */
17+export async function executeRequest(
18+ req: RunRequest,
19+ backend: BackendInfo,
20+ onOutput?: (text: string) => void,
21+ gpu?: GpuProvider,
22+): Promise<RunResult> {
23+ const stages = req.gpu;
24+ if (!stages) return executeBundle(req, backend, onOutput);
25+ const fallback = (note: string): RunResult => ({
26+ ...executeBundle(req, backend, onOutput),
27+ compute: { loop: "numbl", note },
28+ });
29+
30+ if (!gpu) return fallback("WebGPU is not available here, so the time loop ran in numbl");
31+ const shipped = LIBRARY_FILES.find((f) => f.name === stages.loopFile)?.source;
32+ const current = req.bundle.files.find((f) => f.name === stages.loopFile)?.source;
33+ if (shipped !== undefined && current !== shipped)
34+ return fallback(
35+ `${stages.loopFile} was edited, so the time loop ran in numbl, where the edit is live`,
36+ );
37+ let handle: GpuHandle;
38+ try {
39+ handle = await gpu();
40+ } catch (e) {
41+ return fallback(`no WebGPU device (${formatError(e)}), so the time loop ran in numbl`);
42+ }
43+
44+ const t0 = performance.now();
45+ const output: string[] = [];
46+ const say = (text: string) => {
47+ output.push(text);
48+ onOutput?.(text);
49+ };
50+ const libs = req.bundle.files.filter((f) => f.name !== req.bundle.entry);
51+ const fail = (error: string, compute?: ComputeInfo): RunResult => ({
52+ ok: false,
53+ error,
54+ output,
55+ wallMs: performance.now() - t0,
56+ backend,
57+ compute,
58+ });
59+
60+ const setupRun = executeBundle(
61+ {
62+ bundle: {
63+ files: [{ name: "gpu_setup.m", source: stages.setup }, ...libs],
64+ entry: "gpu_setup.m",
65+ resultVariable: "S",
66+ },
67+ },
68+ backend,
69+ say,
70+ );
71+ if (!setupRun.ok) return fail(setupRun.error);
72+ let setup;
73+ try {
74+ setup = fdLoopSetupFromPlain(setupRun.value);
75+ } catch (e) {
76+ return fail(`the setup stage returned an unexpected shape: ${formatError(e)}`);
77+ }
78+
79+ let lastPct = 0;
80+ let rec;
81+ try {
82+ rec = await runFdLoopGpu(handle.device, setup, (done, total) => {
83+ const pct = Math.floor((10 * done) / total);
84+ if (pct > lastPct && done < total) {
85+ lastPct = pct;
86+ say(`gpu: ${done.toLocaleString()} of ${total.toLocaleString()} steps\n`);
87+ }
88+ });
89+ } catch (e) {
90+ return fail(`WebGPU loop: ${formatError(e)}`, { loop: "webgpu", adapter: handle.adapter });
91+ }
92+
93+ const R = fdLoopRecordToPlain(rec, setup, handle.adapter);
94+ const outRun = executeBundle(
95+ {
96+ bundle: {
97+ files: [{ name: "gpu_output.m", source: stages.output }, ...libs],
98+ entry: "gpu_output.m",
99+ resultVariable: req.bundle.resultVariable,
100+ },
101+ inputs: { ...(req.inputs ?? {}), R },
102+ },
103+ backend,
104+ say,
105+ );
106+ const compute: ComputeInfo = {
107+ loop: "webgpu",
108+ adapter: handle.adapter,
109+ gpuMs: rec.gpuMs,
110+ numblMs: setupRun.wallMs + outRun.wallMs,
111+ };
112+ if (!outRun.ok) return fail(outRun.error, compute);
113+ return {
114+ ok: true,
115+ value: outRun.value,
116+ output,
117+ wallMs: performance.now() - t0,
118+ backend,
119+ compute,
120+ };
121+}
packages/core/src/runner/nodeRunner.tsmodified+19−4View file
@@ -1,13 +1,28 @@
11 import type { BackendInfo, RunRequest, RunResult, Runner } from "./types";
2-import { executeBundle, numblVersion } from "./execute";
2+import { numblVersion } from "./execute";
3+import { executeRequest, type GpuProvider } from "./executeGpu";
34
4-/** Runs numbl in the current Node process. Used by tests and scripts; cannot cancel. */
5+/**
6+ * Runs numbl in the current Node process. Used by tests and scripts; cannot
7+ * cancel. Pass `gpu` to run the ported time loop on a WebGPU device (the
8+ * optional `webgpu` package, Google Dawn, provides one under Node).
9+ */
510 export class NodeRunner implements Runner {
6- private readonly info: BackendInfo = { engine: "numbl", version: numblVersion(), optimization: "1", host: "node" };
11+ private readonly info: BackendInfo = {
12+ engine: "numbl",
13+ version: numblVersion(),
14+ optimization: "1",
15+ host: "node",
16+ };
717 readonly backend = Promise.resolve(this.info);
18+ private readonly gpu?: GpuProvider;
19+
20+ constructor(opts: { gpu?: GpuProvider } = {}) {
21+ this.gpu = opts.gpu;
22+ }
823
924 run(req: RunRequest, onOutput?: (text: string) => void): Promise<RunResult> {
10- return Promise.resolve(executeBundle(req, this.info, onOutput));
25+ return executeRequest(req, this.info, onOutput, this.gpu);
1126 }
1227 cancel(): void {
1328 /* synchronous host: nothing in flight between awaits */
packages/core/src/runner/protocol.tsmodified+8−2View file
@@ -1,9 +1,15 @@
11 import type { SourceBundle } from "../generate/bundle";
22 import type { PlainValue } from "./plain";
3-import type { BackendInfo, RunResult } from "./types";
3+import type { BackendInfo, GpuStages, RunResult } from "./types";
44
55 export type ToWorker =
6- | { type: "run"; id: number; bundle: SourceBundle; inputs?: Record<string, PlainValue> }
6+ | {
7+ type: "run";
8+ id: number;
9+ bundle: SourceBundle;
10+ inputs?: Record<string, PlainValue>;
11+ gpu?: GpuStages;
12+ }
713 | { type: "info"; id: number };
814
915 export type FromWorker =
packages/core/src/runner/types.tsmodified+28−0View file
@@ -10,10 +10,36 @@ export interface BackendInfo {
1010 host: "web-worker" | "node";
1111 }
1212
13+/**
14+ * A run whose time loop may leave numbl: the host runs `setup` in numbl, hands
15+ * the struct S to its WebGPU port of `loopFile`, and runs `output` in numbl
16+ * with the loop's record as the workspace variable R. Both scripts carry the
17+ * same params block as the bundle's entry.
18+ */
19+export interface GpuStages {
20+ setup: string;
21+ output: string;
22+ /** The library file whose port runs on the GPU; an edit to it falls back to numbl */
23+ loopFile: string;
24+}
25+
1326 export interface RunRequest {
1427 bundle: SourceBundle;
1528 /** Workspace variables defined before the script runs */
1629 inputs?: Record<string, PlainValue>;
30+ /** Ask for the time loop on the GPU; ignored by a runner without one */
31+ gpu?: GpuStages;
32+}
33+
34+/** Where the time loop of a staged run actually ran, and why, when it was asked to leave numbl. */
35+export interface ComputeInfo {
36+ loop: "webgpu" | "numbl";
37+ adapter?: string;
38+ /** Why the loop stayed in numbl although the GPU was asked for */
39+ note?: string;
40+ /** Time in the GPU loop and in the numbl setup and output stages [ms] */
41+ gpuMs?: number;
42+ numblMs?: number;
1743 }
1844
1945 export interface RunSuccess {
@@ -24,6 +50,7 @@ export interface RunSuccess {
2450 /** Wall time of executeCode, including JIT compilation [ms] */
2551 wallMs: number;
2652 backend: BackendInfo;
53+ compute?: ComputeInfo;
2754 }
2855
2956 export interface RunFailure {
@@ -33,6 +60,7 @@ export interface RunFailure {
3360 output: string[];
3461 wallMs: number;
3562 backend: BackendInfo;
63+ compute?: ComputeInfo;
3664 /** True when the run was cancelled by the host rather than failing */
3765 cancelled?: boolean;
3866 }
packages/core/src/runner/worker.tsmodified+38−8View file
@@ -1,26 +1,56 @@
11 /// <reference lib="webworker" />
22 /**
3- * Web Worker that owns a numbl instance. One run at a time; the host cancels a
4- * run by terminating the worker, which is why nothing here needs to be
5- * interruptible.
3+ * Web Worker that owns a numbl instance and, when the browser exposes WebGPU
4+ * to workers, a GPU device for the ported time loop. One run at a time; the
5+ * host cancels a run by terminating the worker, which is why nothing here
6+ * needs to be interruptible.
67 */
7-import { executeBundle, numblVersion } from "./execute";
8+import { numblVersion } from "./execute";
9+import { executeRequest, type GpuProvider } from "./executeGpu";
10+import { requestGpuDevice, type GpuHandle } from "../gpu/fdLoop";
811 import { collectBuffers } from "./plain";
912 import type { FromWorker, ToWorker } from "./protocol";
1013 import type { BackendInfo } from "./types";
1114
12-const backend: BackendInfo = { engine: "numbl", version: numblVersion(), optimization: "1", host: "web-worker" };
15+const backend: BackendInfo = {
16+ engine: "numbl",
17+ version: numblVersion(),
18+ optimization: "1",
19+ host: "web-worker",
20+};
1321
1422 const post = (msg: FromWorker, transfer: ArrayBuffer[] = []) => self.postMessage(msg, transfer);
1523
16-self.onmessage = (e: MessageEvent<ToWorker>) => {
24+// The device is requested once and kept; a lost device is requested again next time.
25+let gpuHandle: Promise<GpuHandle> | null = null;
26+const gpuProvider: GpuProvider = () => {
27+ if (!gpuHandle) {
28+ gpuHandle = requestGpuDevice().then((h) => {
29+ void h.device.lost.then(() => {
30+ gpuHandle = null;
31+ });
32+ return h;
33+ });
34+ gpuHandle.catch(() => {
35+ gpuHandle = null;
36+ });
37+ }
38+ return gpuHandle;
39+};
40+const hasWebGpu =
41+ typeof navigator !== "undefined" && "gpu" in navigator && !!(navigator as { gpu?: unknown }).gpu;
42+
43+self.onmessage = async (e: MessageEvent<ToWorker>) => {
1744 const msg = e.data;
1845 if (msg.type === "info") {
1946 post({ type: "info", id: msg.id, backend });
2047 return;
2148 }
22- const result = executeBundle({ bundle: msg.bundle, inputs: msg.inputs }, backend, (text) =>
23- post({ type: "output", id: msg.id, text }),
49+ const result = await executeRequest(
50+ { bundle: msg.bundle, inputs: msg.inputs, gpu: msg.gpu },
51+ backend,
52+ (text) => post({ type: "output", id: msg.id, text }),
53+ hasWebGpu ? gpuProvider : undefined,
2454 );
2555 const transfer = result.ok ? collectBuffers(result.value) : [];
2656 post({ type: "done", id: msg.id, result }, transfer);
packages/core/src/runner/workerRunner.tsmodified+1−1View file
@@ -76,7 +76,7 @@ export class WorkerRunner implements Runner {
7676 return new Promise((resolve) => {
7777 const id = this.nextId++;
7878 this.pending.set(id, { resolve, onOutput, output: [], t0: performance.now() });
79- this.ensureWorker().postMessage({ type: "run", id, bundle: req.bundle, inputs: req.inputs } satisfies ToWorker);
79+ this.ensureWorker().postMessage({ type: "run", id, bundle: req.bundle, inputs: req.inputs, gpu: req.gpu } satisfies ToWorker);
8080 });
8181 }
8282
packages/core/src/spec/presets.tsmodified+23−10View file
@@ -64,7 +64,8 @@ export const MATERIAL_PRESETS: MaterialPreset[] = [
6464 id: "1095-spring-steel",
6565 material: { name: "1095 spring steel (bar stock)", E: 205e9, nu: 0.29, rho: 7850 },
6666 damping: { sigma0: 0.5, sigma1: 0.003 },
67- source: "High-carbon steel used for glockenspiel and chime bars; E and rho typical of carbon steel.",
67+ source:
68+ "High-carbon steel used for glockenspiel and chime bars; E and rho typical of carbon steel.",
6869 },
6970 ];
7071
@@ -82,31 +83,43 @@ export const PLATE_PRESETS: PlatePreset[] = [
8283 {
8384 id: "sheet-300x200x2",
8485 name: "300 x 200 x 2 mm sheet",
85- plate: { Lx: 0.3, Ly: 0.2, h: 0.002 },
86+ plate: { shape: "rectangle", Lx: 0.3, Ly: 0.2, h: 0.002 },
8687 note: "A common cut-to-size sheet; the default problem.",
8788 },
8889 {
8990 id: "sheet-12x12x1-8in",
90- name: '12 x 12 x 1/8 in plate',
91- plate: { Lx: 12 * inch, Ly: 12 * inch, h: inch / 8 },
91+ name: "12 x 12 x 1/8 in plate",
92+ plate: { shape: "rectangle", Lx: 12 * inch, Ly: 12 * inch, h: inch / 8 },
9293 note: "Square stock in imperial sizes.",
9394 },
9495 {
9596 id: "sheet-500x300x1",
9697 name: "500 x 300 x 1 mm sheet",
97- plate: { Lx: 0.5, Ly: 0.3, h: 0.001 },
98+ plate: { shape: "rectangle", Lx: 0.5, Ly: 0.3, h: 0.001 },
9899 note: "Thin and large: many modes below 20 kHz, expensive for finite differences.",
99100 },
100101 {
101- id: "cymbal-flat-400x400x1p2",
102- name: "Flat cymbal stand-in, 400 x 400 x 1.2 mm",
103- plate: { Lx: 0.4, Ly: 0.4, h: 0.0012 },
104- note: "A real cymbal is circular, curved and free at the edge; this rectangular, supported plate is only a stand-in for its thickness and span.",
102+ id: "disk-300x2",
103+ name: "Disk, 300 mm diameter x 2 mm",
104+ plate: { shape: "ellipse", Lx: 0.3, Ly: 0.3, h: 0.002 },
105+ note: "A round blank. Free edges and an off-centre strike give the classic Chladni disk; only the finite-difference methods handle the shape.",
106+ },
107+ {
108+ id: "ellipse-300x200x2",
109+ name: "Elliptical plate, 300 x 200 x 2 mm",
110+ plate: { shape: "ellipse", Lx: 0.3, Ly: 0.2, h: 0.002 },
111+ note: "The default sheet cut to an ellipse: same box, same strike, less plate. Compare the two.",
112+ },
113+ {
114+ id: "cymbal-flat-400x1p2",
115+ name: "Flat cymbal stand-in, 400 mm diameter x 1.2 mm",
116+ plate: { shape: "ellipse", Lx: 0.4, Ly: 0.4, h: 0.0012 },
117+ note: "A cymbal is a disk, free at the rim and held at the centre: choose the pinned-centre edge condition. What this flat disk leaves out is the curvature (a real cymbal is a shallow cone with a raised bell) and the taper in thickness.",
105118 },
106119 {
107120 id: "glockenspiel-bar",
108121 name: "Glockenspiel bar, 200 x 32 x 6 mm",
109- plate: { Lx: 0.2, Ly: 0.032, h: 0.006 },
122+ plate: { shape: "rectangle", Lx: 0.2, Ly: 0.032, h: 0.006 },
110123 note: "Real bars are free at both ends and rest on nodal supports; supported or clamped edges give different frequencies. Thickness is a fifth of the width, at the edge of thin-plate validity.",
111124 },
112125 ];
packages/core/src/spec/types.tsmodified+22−3View file
@@ -8,17 +8,36 @@
88 */
99
1010 /** The same condition on all four edges; "free" is the free-free plate of Chladni figures and modal testing. */
11-export type BoundaryCondition = "simply-supported" | "clamped" | "free";
11+/** "pinned-center" is free edges with the node nearest the centre held at w = 0: a cymbal on its stand. */
12+export type BoundaryCondition = "simply-supported" | "clamped" | "free" | "pinned-center";
13+
14+/**
15+ * The outline of the plate within its bounding box [0, Lx] x [0, Ly]: the box
16+ * itself, or the inscribed ellipse (a disk when Lx = Ly). The finite-difference
17+ * methods embed the ellipse in the same grid as the rectangle; the modal method
18+ * has closed-form modes only for the rectangle.
19+ */
20+export type PlateShape = "rectangle" | "ellipse";
1221
1322 export interface PlateGeometry {
14- /** Side along x [m] */
23+ shape: PlateShape;
24+ /** Extent along x [m]: the side of a rectangle, the diameter of an ellipse */
1525 Lx: number;
16- /** Side along y [m] */
26+ /** Extent along y [m] */
1727 Ly: number;
1828 /** Thickness [m] */
1929 h: number;
2030 }
2131
32+/** Whether a point of the bounding box lies on the plate. */
33+export function onPlate(plate: PlateGeometry, x: number, y: number): boolean {
34+ if (!(x >= 0 && x <= plate.Lx && y >= 0 && y <= plate.Ly)) return false;
35+ if (plate.shape !== "ellipse") return true;
36+ const a = plate.Lx / 2;
37+ const b = plate.Ly / 2;
38+ return ((x - a) / a) ** 2 + ((y - b) / b) ** 2 <= 1;
39+}
40+
2241 export interface Material {
2342 name: string;
2443 /** Young's modulus [Pa] */
packages/core/src/spec/validate.tsmodified+19−9View file
@@ -1,4 +1,4 @@
1-import type { ProblemSpec } from "./types";
1+import { onPlate, type ProblemSpec } from "./types";
22
33 export interface ValidationIssue {
44 /** Dotted path into the spec, e.g. "plate.h" */
@@ -28,9 +28,9 @@ export function validateSpec(spec: ProblemSpec): ValidationIssue[] {
2828 if (!(material.rho > 0)) err("material.rho", "must be positive");
2929 if (damping.sigma0 < 0 || damping.sigma1 < 0) err("damping", "loss rates cannot be negative");
3030
31- const inside = (x: number, y: number) => x >= 0 && x <= plate.Lx && y >= 0 && y <= plate.Ly;
32- if (!inside(excitation.x, excitation.y)) err("excitation", "strike point lies outside the plate");
33- if (!inside(probe.x, probe.y)) err("probe", "probe lies outside the plate");
31+ const where = plate.shape === "ellipse" ? "outside the ellipse" : "outside the plate";
32+ if (!onPlate(plate, excitation.x, excitation.y)) err("excitation", `strike point lies ${where}`);
33+ if (!onPlate(plate, probe.x, probe.y)) err("probe", `probe lies ${where}`);
3434 if (!(excitation.radius > 0)) err("excitation.radius", "footprint radius must be positive");
3535 if (excitation.type === "impulse" && !(excitation.duration > 0))
3636 err("excitation.duration", "contact time must be positive");
@@ -43,14 +43,24 @@ export function validateSpec(spec: ProblemSpec): ValidationIssue[] {
4343 if (modal.fmax > output.fs / 2)
4444 warn("modal.fmax", "modes above half the sample rate alias in the probe series and audio");
4545 if (fd.nx < 5 || fd.ny < 5) err("fd", "finite-difference grid needs at least 5 nodes per side");
46- if (fd.cfl > 1) warn("fd.cfl", "above the stability limit: the finite-difference run will blow up");
47- if (implicit.nx < 5 || implicit.ny < 5) err("implicit", "implicit grid needs at least 5 nodes per side");
46+ if (fd.cfl > 1)
47+ warn("fd.cfl", "above the stability limit: the finite-difference run will blow up");
48+ if (implicit.nx < 5 || implicit.ny < 5)
49+ err("implicit", "implicit grid needs at least 5 nodes per side");
4850 if (!(implicit.steps_per_sample > 0)) err("implicit.steps_per_sample", "must be positive");
49- if (implicit.theta < 0.25) warn("implicit.theta", "below 1/4 the implicit scheme is only conditionally stable");
50- if (implicit.steps_per_sample < 1) warn("implicit.steps_per_sample", "fewer than one step per sample detunes modes near half the sample rate by tens of percent");
51+ if (implicit.theta < 0.25)
52+ warn("implicit.theta", "below 1/4 the implicit scheme is only conditionally stable");
53+ if (implicit.steps_per_sample < 1)
54+ warn(
55+ "implicit.steps_per_sample",
56+ "fewer than one step per sample detunes modes near half the sample rate by tens of percent",
57+ );
5158 const dxFd = plate.Lx / (fd.nx - 1);
5259 const dyFd = plate.Ly / (fd.ny - 1);
5360 if (excitation.radius < Math.max(dxFd, dyFd))
54- warn("excitation.radius", "footprint smaller than a finite-difference cell: the load lands on one node");
61+ warn(
62+ "excitation.radius",
63+ "footprint smaller than a finite-difference cell: the load lands on one node",
64+ );
5565 return out;
5666 }
packages/core/test/ellipse.test.tsadded+262−0View file
@@ -0,0 +1,262 @@
1+import { describe, expect, it } from "vitest";
2+import {
3+ LIBRARY_FILES,
4+ buildCompareBundle,
5+ isPlainStruct,
6+ methodBlockers,
7+ getMethod,
8+ parseCompareResult,
9+ validateSpec,
10+} from "../src/index";
11+import { maxAbs, runner, solve, testSpec } from "./helpers";
12+
13+/** Run a MATLAB snippet against the library and return the struct it leaves in `out`. */
14+async function script(source: string) {
15+ const run = await runner.run({
16+ bundle: {
17+ files: [{ name: "c.m", source }, ...LIBRARY_FILES],
18+ entry: "c.m",
19+ resultVariable: "out",
20+ },
21+ });
22+ if (!run.ok) throw new Error(`${run.error}\n${run.output.join("")}`);
23+ if (!isPlainStruct(run.value)) throw new Error("unexpected shape");
24+ return run.value.fields as Record<string, number>;
25+}
26+
27+function ellipseSpec() {
28+ const s = testSpec();
29+ s.plate.shape = "ellipse";
30+ return s;
31+}
32+
33+describe("the plate's footprint on the grid", () => {
34+ it("is the whole grid for a rectangle and a staircase for an ellipse, with every node carrying area", async () => {
35+ const out = await script(`
36+ r = plate_domain(struct('shape', 'rectangle', 'Lx', 0.3, 'Ly', 0.2), 19, 13);
37+ out.r_inside = nnz(r.inside); out.r_interior = nnz(r.interior); out.r_boundary = nnz(r.boundary);
38+ out.r_m_corner = r.m(1, 1); out.r_m_edge = r.m(1, 7); out.r_m_in = r.m(5, 7); out.r_full3 = nnz(r.full3);
39+ e = plate_domain(struct('shape', 'ellipse', 'Lx', 0.3, 'Ly', 0.2), 31, 21);
40+ out.e_inside = nnz(e.inside); out.e_area = nnz(e.inside) / (31 * 21);
41+ out.e_min_m = min(e.m(e.inside)); out.e_max_m = max(e.m(:));
42+ out.e_corner = e.inside(1, 1); out.e_centre = e.inside(11, 16);
43+ out.e_near_ok = all(e.full3(e.near(e.inside)));
44+ out.e_kappa_min = min(e.kappa(e.boundary)); out.e_kappa_max = max(e.kappa(e.boundary));
45+ nrm = sqrt(e.normal_x(e.boundary).^2 + e.normal_y(e.boundary).^2);
46+ out.e_normal_unit = max(abs(nrm - 1));
47+ % a sliver: the disk on a coarse grid still has full stencils
48+ d = plate_domain(struct('shape', 'ellipse', 'Lx', 1, 'Ly', 1), 9, 9);
49+ out.d_full3 = nnz(d.full3);
50+ `);
51+ expect(out.r_inside).toBe(19 * 13);
52+ expect(out.r_interior).toBe(17 * 11);
53+ expect(out.r_boundary).toBe(19 * 13 - 17 * 11);
54+ expect(out.r_m_corner).toBe(0.25);
55+ expect(out.r_m_edge).toBe(0.5);
56+ expect(out.r_m_in).toBe(1);
57+ expect(out.r_full3).toBe(17 * 11);
58+ // pi/4 of the box, to within the staircase
59+ expect(Math.abs(out.e_area - Math.PI / 4)).toBeLessThan(0.05);
60+ expect(out.e_min_m).toBeGreaterThanOrEqual(0.25);
61+ expect(out.e_max_m).toBe(1);
62+ expect(out.e_corner).toBe(false);
63+ expect(out.e_centre).toBe(true);
64+ expect(out.e_near_ok).toBe(true);
65+ // curvature of a 300 x 200 ellipse lies between b/a^2 and a/b^2
66+ expect(out.e_kappa_min).toBeGreaterThan(0.1 / 0.15 ** 2 - 1e-9);
67+ expect(out.e_kappa_max).toBeLessThan(0.15 / 0.1 ** 2 + 1e-9);
68+ expect(out.e_normal_unit).toBeLessThan(1e-12);
69+ expect(out.d_full3).toBeGreaterThan(0);
70+ });
71+
72+ it("keeps the strike and the probe on the plate", () => {
73+ const spec = ellipseSpec();
74+ expect(validateSpec(spec).some((i) => i.severity === "error")).toBe(false);
75+ // inside the box, outside the ellipse
76+ spec.excitation.x = 0.01;
77+ spec.excitation.y = 0.01;
78+ const issues = validateSpec(spec);
79+ expect(issues.some((i) => i.path === "excitation" && /ellipse/.test(i.message))).toBe(true);
80+ });
81+});
82+
83+describe("circular plate frequencies", () => {
84+ // omega a^2 sqrt(rho h / D) for the disk of radius a, nu = 0.3, from the Bessel-function
85+ // frequency equations (clamped: J_n I_{n+1} + I_n J_{n+1} = 0; simply supported:
86+ // J_{n+1}/J_n + I_{n+1}/I_n = 2 lambda / (1 - nu)) and, for free edges, the tabulated values
87+ // (Leissa, Vibration of Plates, 1969). The staircase edge is first order, so a 25-node grid
88+ // is expected within a few percent.
89+ it("clamped: 10.216", async () => {
90+ const out = await script(`
91+ dom = plate_domain(struct('shape', 'ellipse', 'Lx', 2, 'Ly', 2), 25, 25);
92+ [~, K, M] = plate_matrices(dom, 'clamped', 0.3);
93+ ev = sort(eig(full(K)));
94+ out.f1 = sqrt(ev(1)); out.f2 = sqrt(ev(2)); out.f3 = sqrt(ev(3));
95+ out.ratio = ev(end) / (4/dom.dx^2 + 4/dom.dy^2)^2;
96+ out.n = size(K, 1);
97+ `);
98+ expect(Math.abs(out.f1 - 10.216) / 10.216).toBeLessThan(0.02);
99+ // the (1, 0) pair, degenerate on the true disk
100+ expect(Math.abs(out.f2 - 21.26) / 21.26).toBeLessThan(0.03);
101+ expect(Math.abs(out.f3 - 21.26) / 21.26).toBeLessThan(0.03);
102+ // the explicit step bound derived for the interior stencil still holds
103+ expect(out.ratio).toBeLessThan(1);
104+ expect(out.ratio).toBeGreaterThan(0.9);
105+ });
106+
107+ it("simply supported: 4.935 with the curvature term, the paradox value 5.783 without it", async () => {
108+ const out = await script(`
109+ dom = plate_domain(struct('shape', 'ellipse', 'Lx', 2, 'Ly', 2), 25, 25);
110+ [L, K] = plate_matrices(dom, 'simply-supported', 0.3);
111+ ev = sort(eig(full(K)));
112+ out.f1 = sqrt(ev(1)); out.f2 = sqrt(ev(2));
113+ out.minev = ev(1) / ev(end);
114+ ev0 = sort(eig(full(L * L)));
115+ out.g1 = sqrt(ev0(1));
116+ out.symm = norm(full(K - K'), 1) / norm(full(K), 1);
117+ `);
118+ expect(out.symm).toBeLessThan(1e-12);
119+ expect(out.minev).toBeGreaterThan(0);
120+ expect(Math.abs(out.f1 - 4.935) / 4.935).toBeLessThan(0.03);
121+ expect(Math.abs(out.f2 - 13.898) / 13.898).toBeLessThan(0.03);
122+ // L*L alone converges to the plate with zero Laplacian on the edge, J_0's first zero squared
123+ expect(Math.abs(out.g1 - 5.783) / 5.783).toBeLessThan(0.03);
124+ expect(out.g1 / out.f1).toBeGreaterThan(1.12);
125+ });
126+
127+ it("free: 5.358 and 9.003, with exactly three rigid-body modes", async () => {
128+ const out = await script(`
129+ dom = plate_domain(struct('shape', 'ellipse', 'Lx', 2, 'Ly', 2), 25, 25);
130+ [L, K, M] = plate_matrices(dom, 'free', 0.3);
131+ Dh = diag(1 ./ sqrt(full(diag(M)))); S = Dh * full(K) * Dh; S = (S + S') / 2;
132+ ev = sort(eig(S));
133+ out.zero3 = max(abs(ev(1:3))) / ev(end);
134+ out.fourth = ev(4) / ev(end);
135+ out.f1 = sqrt(ev(4)); out.f2 = sqrt(ev(5)); out.f3 = sqrt(ev(6));
136+ out.ratio = ev(end) / (4/dom.dx^2 + 4/dom.dy^2)^2;
137+ X = dom.X(dom.inside); Y = dom.Y(dom.inside);
138+ out.nullx = max(abs(K * (X - 1))) * dom.dx^4; out.nully = max(abs(K * (Y - 1))) * dom.dx^4;
139+ out.lnull = max(abs(L * ones(size(X))));
140+ `);
141+ expect(out.zero3).toBeLessThan(1e-9);
142+ expect(out.fourth).toBeGreaterThan(1e-6);
143+ expect(out.nullx).toBeLessThan(1e-9);
144+ expect(out.nully).toBeLessThan(1e-9);
145+ expect(out.lnull).toBeLessThan(1e-9);
146+ // the (2, 0) pair is split by the staircase's anisotropy; its mean is the disk's value
147+ expect(Math.abs((out.f1 + out.f2) / 2 - 5.358) / 5.358).toBeLessThan(0.03);
148+ expect(Math.abs(out.f3 - 9.003) / 9.003).toBeLessThan(0.03);
149+ expect(out.ratio).toBeLessThan(1);
150+ });
151+});
152+
153+describe("solving on an ellipse", () => {
154+ it.each(["simply-supported", "clamped", "free", "pinned-center"] as const)(
155+ "runs both finite-difference schemes with %s edges and they agree",
156+ async (bc) => {
157+ const spec = ellipseSpec();
158+ spec.bc = bc;
159+ spec.implicit = {
160+ nx: spec.fd.nx,
161+ ny: spec.fd.ny,
162+ duration: spec.fd.duration,
163+ steps_per_sample: 2,
164+ theta: 0.25,
165+ };
166+ const exp = await solve(spec, "fd");
167+ const imp = await solve(spec, "implicit");
168+ for (const { result } of [exp, imp]) {
169+ expect(result.info.stable).toBe(true);
170+ expect(result.info.extra.shape).toBe("ellipse");
171+ expect(result.probe.w.every(Number.isFinite)).toBe(true);
172+ expect(maxAbs(result.probe.w)).toBeGreaterThan(0);
173+ expect(result.info.notes.some((n) => /staircase/.test(n))).toBe(true);
174+ // outside the ellipse the fields are NaN, and the same nodes in every frame and field
175+ const { nx, ny, nf, data } = result.frames.w;
176+ const nanPer = Array.from(
177+ { length: nf },
178+ (_, k) =>
179+ data.subarray(k * nx * ny, (k + 1) * nx * ny).filter((v) => Number.isNaN(v)).length,
180+ );
181+ expect(nanPer[0]).toBeGreaterThan(0.1 * nx * ny);
182+ expect(nanPer[0]).toBeLessThan(0.3 * nx * ny);
183+ expect(new Set(nanPer).size).toBe(1);
184+ expect(result.frames.vm.data.filter((v) => Number.isNaN(v)).length).toBe(nanPer[0] * nf);
185+ expect(result.envelope.vm.data.filter((v) => Number.isNaN(v)).length).toBe(nanPer[0]);
186+ // the corners of the box are off the plate, the centre is on it
187+ expect(data[0]).toBeNaN();
188+ expect(data[ny - 1 + (nx - 1) * ny]).toBeNaN();
189+ expect(Number.isFinite(data[Math.floor(ny / 2) + Math.floor(nx / 2) * ny])).toBe(true);
190+ // the probe was kept where it was asked for (0.2, 0.12), on the plate
191+ expect(Math.abs(result.probe.x - spec.probe.x)).toBeLessThan(
192+ 0.5 * (spec.plate.Lx / (nx - 1)),
193+ );
194+ expect(Math.abs(result.probe.y - spec.probe.y)).toBeLessThan(
195+ 0.5 * (spec.plate.Ly / (ny - 1)),
196+ );
197+ }
198+ expect(Number(imp.result.info.extra.nunknowns)).toBeLessThan(spec.fd.nx * spec.fd.ny);
199+ expect(Number(exp.result.info.extra.nunknowns)).toBe(Number(imp.result.info.extra.nunknowns));
200+ const cmp = await runner.run({
201+ bundle: buildCompareBundle(),
202+ inputs: { A: exp.run.ok ? exp.run.value : 0, B: imp.run.ok ? imp.run.value : 0 },
203+ });
204+ expect(cmp.ok).toBe(true);
205+ if (!cmp.ok) return;
206+ const c = parseCompareResult(cmp.value);
207+ expect(c.w_rel_l2).toBeLessThan(0.03);
208+ expect(c.vm_rel_l2).toBeLessThan(0.1);
209+ expect(c.frames_rel_l2.length).toBeGreaterThan(0);
210+ expect(c.frames_rel_l2.every(Number.isFinite)).toBe(true);
211+ },
212+ );
213+
214+ it("is stiffer than the rectangle of the same box, and diverges above the CFL limit like it", async () => {
215+ // A clamped plate on a smaller domain has a higher fundamental (domain monotonicity of the
216+ // Dirichlet problem), so the ellipse inscribed in the box must beat the box itself.
217+ const out = await script(`
218+ r = plate_domain(struct('shape', 'rectangle', 'Lx', 0.3, 'Ly', 0.2), 25, 17);
219+ e = plate_domain(struct('shape', 'ellipse', 'Lx', 0.3, 'Ly', 0.2), 25, 17);
220+ [~, Kr] = plate_matrices(r, 'clamped', 0.3); [~, Ke] = plate_matrices(e, 'clamped', 0.3);
221+ out.rect = min(eig(full(Kr))); out.ell = min(eig(full(Ke)));
222+ `);
223+ expect(out.ell).toBeGreaterThan(out.rect);
224+ // The Rayleigh quotient of the trial function (1 - x^2/a^2 - y^2/b^2)^2 is an upper bound,
225+ // omega^2 rho h / D = (40/3) (3/a^4 + 2/(a^2 b^2) + 3/b^4): 10.33 against 10.216 for the disk,
226+ // 17.39 for a/b = 1.5. The true value lies about a percent below it, and the staircase on a
227+ // 25-node grid a further percent or two below that (the disk gives 10.21 at this resolution).
228+ const a = 0.15;
229+ const b = 0.1;
230+ const bound = Math.sqrt((40 / 3) * (3 / a ** 4 + 2 / (a * b) ** 2 + 3 / b ** 4)) * a ** 2;
231+ expect(bound).toBeCloseTo(17.39, 1);
232+ const lam = Math.sqrt(out.ell) * a ** 2;
233+ expect(lam).toBeLessThan(bound);
234+ expect(lam).toBeGreaterThan(0.94 * bound);
235+ const ell = ellipseSpec();
236+ ell.fd.cfl = 1.04;
237+ const blown = await solve(ell, "fd");
238+ expect(blown.result.info.stable).toBe(false);
239+ expect(blown.result.probe.w.every(Number.isFinite)).toBe(true);
240+ });
241+
242+ it("is refused by the modal method, in the registry and in the solver", async () => {
243+ const spec = ellipseSpec();
244+ expect(methodBlockers(getMethod("modal"), spec).some((b) => /ellipse/.test(b))).toBe(true);
245+ expect(methodBlockers(getMethod("fd"), spec)).toEqual([]);
246+ expect(methodBlockers(getMethod("implicit"), spec)).toEqual([]);
247+ await expect(solve(spec, "modal")).rejects.toThrow(/rectangle/);
248+ });
249+
250+ it("puts the load and the probe on the plate even when the strike is near the edge", async () => {
251+ const spec = ellipseSpec();
252+ // just inside the ellipse near its right end, footprint spilling over the edge
253+ spec.excitation.x = 0.28;
254+ spec.excitation.y = 0.1;
255+ spec.probe.x = 0.27;
256+ spec.probe.y = 0.1;
257+ const { result } = await solve(spec, "fd");
258+ expect(result.info.stable).toBe(true);
259+ expect(maxAbs(result.probe.w)).toBeGreaterThan(0);
260+ expect(Number.isFinite(result.probe.x)).toBe(true);
261+ });
262+});
packages/core/test/free.test.tsmodified+114−5View file
@@ -20,7 +20,8 @@ describe("free edges", () => {
2020 it("assembles a symmetric stiffness whose null space is exactly the rigid-body motions and whose interior is the 13-point stencil", async () => {
2121 const out = await script(`
2222 nx = 19; ny = 13; dx = 0.3/(nx-1); dy = 0.2/(ny-1); nu = 0.3;
23- [L, K, M, idx, mx, my] = plate_matrices(nx, ny, dx, dy, 'free', nu);
23+ dom = plate_domain(struct('shape', 'rectangle', 'Lx', 0.3, 'Ly', 0.2), nx, ny);
24+ [L, K, M, idx] = plate_matrices(dom, 'free', nu);
2425 x = linspace(0, 0.3, nx); y = linspace(0, 0.2, ny); [X, Y] = meshgrid(x, y);
2526 out.symm = norm(full(K - K'), 1) / norm(full(K), 1);
2627 out.null1 = max(abs(K * ones(nx*ny, 1)));
@@ -36,7 +37,7 @@ describe("free edges", () => {
3637 ref = Bw(2:end-1, 2:end-1);
3738 out.interior = norm(inner(:) - ref(:)) / norm(ref(:));
3839 out.m_inner = Mv(5, 7); out.m_edge = Mv(1, 7); out.m_corner = Mv(1, 1);
39- out.nunk = numel(idx); out.mx = mx; out.my = my;
40+ out.nunk = numel(idx);
4041 `);
4142 expect(out.symm).toBeLessThan(1e-12);
4243 expect(out.null1).toBeLessThan(1e-6);
@@ -48,8 +49,6 @@ describe("free edges", () => {
4849 expect(out.m_edge).toBe(0.5);
4950 expect(out.m_corner).toBe(0.25);
5051 expect(out.nunk).toBe(19 * 13);
51- expect(out.mx).toBe(19);
52- expect(out.my).toBe(13);
5352 });
5453
5554 it("reproduces the frequencies of a completely free square plate (Leissa, nu = 0.3)", async () => {
@@ -59,7 +58,8 @@ describe("free edges", () => {
5958 // should land within a couple of percent, and the three rigid modes must sit at zero.
6059 const out = await script(`
6160 n = 17; a = 1; h = a/(n-1); nu = 0.3;
62- [~, K, M] = plate_matrices(n, n, h, h, 'free', nu);
61+ dom = plate_domain(struct('shape', 'rectangle', 'Lx', a, 'Ly', a), n, n);
62+ [~, K, M] = plate_matrices(dom, 'free', nu);
6363 Dh = diag(1 ./ sqrt(full(diag(M)))); S = Dh * full(K) * Dh;
6464 S = (S + S') / 2;
6565 ev = sort(eig(S));
@@ -166,3 +166,112 @@ describe("free edges", () => {
166166 await expect(solve(spec, "modal")).rejects.toThrow(/simply supported/);
167167 });
168168 });
169+
170+describe("pinned centre", () => {
171+ it("removes the centre node from the system and leaves exactly the two rotations about it at zero frequency", async () => {
172+ const out = await script(`
173+ n = 15; a = 1; h = a/(n-1); nu = 0.3;
174+ dom = plate_domain(struct('shape', 'rectangle', 'Lx', a, 'Ly', a), n, n);
175+ [L, K, M, idx, fixed] = plate_matrices(dom, 'pinned-center', nu);
176+ ic = (floor(n/2) + 1) + floor(n/2) * n;
177+ out.fixed = fixed; out.ic = ic;
178+ out.row = max(abs(K(ic, :))); out.col = max(abs(K(:, ic))); out.lrow = max(abs(L(ic, :)));
179+ out.mdiag = full(M(ic, ic));
180+ x = linspace(0, a, n); [X, Y] = meshgrid(x, x);
181+ out.rotx = max(abs(K * (X(:) - X(ic)))) * h^4; out.roty = max(abs(K * (Y(:) - Y(ic)))) * h^4;
182+ out.const = max(abs(K * ones(n*n, 1))) * h^4;
183+ % the pinned node's own row and column are zero, so drop it before looking at the spectrum
184+ keep = true(n*n, 1); keep(ic) = false;
185+ Dh = diag(1 ./ sqrt(full(diag(M)))); S = Dh * full(K) * Dh; S = (S + S') / 2; S = S(keep, keep);
186+ ev = sort(eig(S));
187+ out.zero2 = max(abs(ev(1:2))) / ev(end);
188+ out.third = ev(3) / ev(end);
189+ out.f1 = sqrt(ev(3)) * a^2; out.f2 = sqrt(ev(4)) * a^2;
190+ [~, Kf, Mf] = plate_matrices(dom, 'free', nu);
191+ Dhf = diag(1 ./ sqrt(full(diag(Mf)))); Sf = Dhf * full(Kf) * Dhf; Sf = (Sf + Sf') / 2;
192+ evf = sort(eig(Sf));
193+ out.free_f1 = sqrt(evf(4)) * a^2;
194+ `);
195+ expect(out.fixed).toBe(out.ic);
196+ expect(out.row).toBe(0);
197+ expect(out.col).toBe(0);
198+ expect(out.lrow).toBe(0);
199+ expect(out.mdiag).toBe(1);
200+ expect(out.rotx).toBeLessThan(1e-9);
201+ expect(out.roty).toBeLessThan(1e-9);
202+ expect(out.const).toBeGreaterThan(1e-3);
203+ expect(out.zero2).toBeLessThan(1e-9);
204+ // the pin removes the translation, so the third eigenvalue is an elastic mode, not a zero
205+ expect(out.third).toBeGreaterThan(1e-6);
206+ // Rayleigh's interlacing theorem: one constraint lifts each eigenvalue to at most the next one,
207+ // so the first elastic mode of the pinned plate (the rim moving in phase, centre held) lies
208+ // between the free plate's third zero and its fundamental
209+ expect(out.f1).toBeGreaterThan(0);
210+ expect(out.f1).toBeLessThan(out.free_f1);
211+ // the free plate's first elastic mode twists about both diagonals and is zero at the centre
212+ // node, so it is an eigenvector of the pinned plate as well, at the same frequency
213+ expect(Math.abs(out.f2 - out.free_f1) / out.free_f1).toBeLessThan(1e-8);
214+ });
215+
216+ it("holds the centre at zero through both schemes while the edges move, and reports the support", async () => {
217+ const spec = testSpec();
218+ spec.bc = "pinned-center";
219+ spec.implicit = {
220+ nx: spec.fd.nx,
221+ ny: spec.fd.ny,
222+ duration: spec.fd.duration,
223+ steps_per_sample: 2,
224+ theta: 0.25,
225+ };
226+ for (const id of ["fd", "implicit"]) {
227+ const { result, run } = await solve(spec, id);
228+ expect(result.info.stable).toBe(true);
229+ expect(result.info.notes.some((n) => /pinned/.test(n))).toBe(true);
230+ const { nx, ny, nf, data } = result.frames.w;
231+ const jc = Math.floor(nx / 2);
232+ const ic = Math.floor(ny / 2);
233+ expect(result.info.extra.support_x).toBeCloseTo(result.grid.x[jc], 12);
234+ expect(result.info.extra.support_y).toBeCloseTo(result.grid.y[ic], 12);
235+ let peak = 0;
236+ let edgePeak = 0;
237+ for (let k = 0; k < nf; k++) {
238+ expect(data[ic + jc * ny + k * nx * ny]).toBe(0);
239+ for (let j = 0; j < nx; j++) {
240+ edgePeak = Math.max(edgePeak, Math.abs(data[0 + j * ny + k * nx * ny]));
241+ for (let i = 0; i < ny; i++)
242+ peak = Math.max(peak, Math.abs(data[i + j * ny + k * nx * ny]));
243+ }
244+ }
245+ expect(peak).toBeGreaterThan(0);
246+ expect(edgePeak).toBeGreaterThan(0.05 * peak);
247+ if (id === "fd") (spec as unknown as { _fd: typeof run })._fd = run;
248+ }
249+ });
250+
251+ it("agrees between the implicit and explicit schemes", async () => {
252+ const spec = testSpec();
253+ spec.bc = "pinned-center";
254+ spec.implicit = {
255+ nx: spec.fd.nx,
256+ ny: spec.fd.ny,
257+ duration: spec.fd.duration,
258+ steps_per_sample: 2,
259+ theta: 0.25,
260+ };
261+ const imp = await solve(spec, "implicit");
262+ const exp = await solve(spec, "fd");
263+ const cmp = await runner.run({
264+ bundle: buildCompareBundle(),
265+ inputs: { A: exp.run.ok ? exp.run.value : 0, B: imp.run.ok ? imp.run.value : 0 },
266+ });
267+ expect(cmp.ok).toBe(true);
268+ if (!cmp.ok) return;
269+ expect(parseCompareResult(cmp.value).w_rel_l2).toBeLessThan(0.03);
270+ });
271+
272+ it("is refused by the modal method", async () => {
273+ const spec = testSpec();
274+ spec.bc = "pinned-center";
275+ await expect(solve(spec, "modal")).rejects.toThrow(/simply supported/);
276+ });
277+});
packages/core/test/gpu.test.tsadded+177−0View file
@@ -0,0 +1,177 @@
1+import { describe, expect, it } from "vitest";
2+import { NodeRunner } from "../src/node";
3+import {
4+ FD_METHOD,
5+ LIBRARY_FILES,
6+ buildGpuStages,
7+ buildRunBundle,
8+ parsePlateResult,
9+ withEdits,
10+ type ProblemSpec,
11+ type RunResult,
12+} from "../src/index";
13+import { requestGpuDevice, type GpuHandle } from "../src/gpu/fdLoop";
14+import { maxAbs, relL2, testSpec } from "./helpers";
15+
16+/**
17+ * Desktop WebGPU for Node through the optional `webgpu` package (Google Dawn):
18+ * `npm i --no-save webgpu`. Without the package, or on a machine without a GPU
19+ * adapter, the suite is skipped rather than failed, so CI can run it
20+ * opportunistically.
21+ */
22+async function dawn(): Promise<GpuHandle | null> {
23+ try {
24+ const specifier = "webgpu";
25+ const mod = (await import(specifier)) as {
26+ create: (flags: string[]) => GPU;
27+ globals: Record<string, unknown>;
28+ };
29+ Object.assign(globalThis, mod.globals);
30+ Object.defineProperty(globalThis, "navigator", {
31+ value: { gpu: mod.create([]) },
32+ configurable: true,
33+ writable: true,
34+ });
35+ return await requestGpuDevice();
36+ } catch (e) {
37+ console.log(
38+ `gpu.test: no WebGPU device under Node (${e instanceof Error ? e.message : e}); skipping`,
39+ );
40+ return null;
41+ }
42+}
43+
44+const handle = await dawn();
45+const cpu = new NodeRunner();
46+const gpu = new NodeRunner({ gpu: () => Promise.resolve(handle!) });
47+
48+async function both(spec: ProblemSpec, edits: Record<string, string> = {}) {
49+ const stages = buildGpuStages(spec, FD_METHOD)!;
50+ const g = await gpu.run({
51+ bundle: withEdits(buildRunBundle(spec, FD_METHOD, { gpu: true }), edits),
52+ gpu: stages,
53+ });
54+ const c = await cpu.run({ bundle: withEdits(buildRunBundle(spec, FD_METHOD), edits) });
55+ if (!g.ok) throw new Error(`gpu run failed: ${g.error}\n${g.output.join("")}`);
56+ if (!c.ok) throw new Error(`cpu run failed: ${c.error}\n${c.output.join("")}`);
57+ return { g, c, rg: parsePlateResult(g.value), rc: parsePlateResult(c.value) };
58+}
59+
60+describe.skipIf(!handle)("explicit finite differences on WebGPU", () => {
61+ it.each(["simply-supported", "clamped", "free", "pinned-center"] as const)(
62+ "reproduces the double-precision loop for %s edges",
63+ async (bc) => {
64+ const spec = testSpec();
65+ spec.bc = bc;
66+ const { g, rg, rc } = await both(spec);
67+ expect(g.compute?.loop).toBe("webgpu");
68+ expect(rg.info.extra.loop).toBe("webgpu");
69+ expect(rg.info.stable).toBe(true);
70+ expect(rg.info.notes.some((n) => /GPU/.test(n))).toBe(true);
71+ expect(rg.t.length).toBe(rc.t.length);
72+ expect(rg.frames.t.length).toBe(rc.frames.t.length);
73+ const eW = relL2(rc.probe.w, rg.probe.w);
74+ const eF = relL2(rc.frames.w.data, rg.frames.w.data);
75+ const eVm = relL2(rc.probe.vm, rg.probe.vm);
76+ const eA = relL2(rc.audio.p, rg.audio.p);
77+ const eEnv = relL2(rc.envelope.vm.data, rg.envelope.vm.data);
78+ console.log(
79+ `gpu.test ${bc}: rel L2 probe w ${eW.toExponential(2)}, frames ${eF.toExponential(2)}, probe vm ${eVm.toExponential(2)}, audio ${eA.toExponential(2)}, envelope ${eEnv.toExponential(2)}, ${rg.info.nsteps} steps in ${g.compute?.gpuMs?.toFixed(0)} ms on the GPU`,
80+ );
81+ expect(eW).toBeLessThan(1e-3);
82+ expect(eF).toBeLessThan(1e-3);
83+ expect(eVm).toBeLessThan(1e-2);
84+ expect(eA).toBeLessThan(1e-2);
85+ expect(eEnv).toBeLessThan(1e-2);
86+ expect(maxAbs(rg.probe.w)).toBeGreaterThan(0);
87+ },
88+ );
89+
90+ it.each(["simply-supported", "clamped", "free", "pinned-center"] as const)(
91+ "reproduces the double-precision loop on an ellipse with %s edges",
92+ async (bc) => {
93+ const spec = testSpec();
94+ spec.plate.shape = "ellipse";
95+ spec.bc = bc;
96+ const { g, rg, rc } = await both(spec);
97+ expect(g.compute?.loop).toBe("webgpu");
98+ expect(rg.info.stable).toBe(true);
99+ expect(rg.frames.w.data.filter((v) => Number.isNaN(v)).length).toBe(
100+ rc.frames.w.data.filter((v) => Number.isNaN(v)).length,
101+ );
102+ const fin = (a: Float64Array) => a.filter(Number.isFinite);
103+ const eW = relL2(rc.probe.w, rg.probe.w);
104+ const eF = relL2(fin(rc.frames.w.data), fin(rg.frames.w.data));
105+ const eVm = relL2(rc.probe.vm, rg.probe.vm);
106+ const eA = relL2(rc.audio.p, rg.audio.p);
107+ console.log(
108+ `gpu.test ellipse ${bc}: rel L2 probe w ${eW.toExponential(2)}, frames ${eF.toExponential(2)}, probe vm ${eVm.toExponential(2)}, audio ${eA.toExponential(2)}`,
109+ );
110+ expect(eW).toBeLessThan(1e-3);
111+ expect(eF).toBeLessThan(1e-3);
112+ expect(eVm).toBeLessThan(1e-2);
113+ expect(eA).toBeLessThan(1e-2);
114+ expect(maxAbs(rg.probe.w)).toBeGreaterThan(0);
115+ },
116+ );
117+
118+ it("stops at the blow-up where the double-precision loop does", async () => {
119+ const spec = testSpec();
120+ spec.fd.cfl = 1.03;
121+ const { rg, rc } = await both(spec);
122+ expect(rg.info.stable).toBe(false);
123+ expect(rc.info.stable).toBe(false);
124+ expect(rg.info.diverged_step).toBeGreaterThan(0);
125+ // the unstable mode grows from rounding noise once the strike is over, so the two runs need not agree to the step
126+ expect(Math.abs(rg.info.diverged_step - rc.info.diverged_step)).toBeLessThanOrEqual(
127+ Math.max(10, 0.1 * rc.info.nsteps),
128+ );
129+ expect(rg.probe.w.every(Number.isFinite)).toBe(true);
130+ expect(rg.frames.w.data.every(Number.isFinite)).toBe(true);
131+ expect(rg.frames.t.length).toBe(rc.frames.t.length);
132+ expect(rg.info.notes.some((n) => /diverged/.test(n))).toBe(true);
133+ });
134+
135+ it("falls back to numbl when the loop file is edited, and says so", async () => {
136+ const spec = testSpec();
137+ const loop = LIBRARY_FILES.find((f) => f.name === "plate_fd_loop.m")!.source;
138+ const { g, rg, rc } = await both(spec, {
139+ "plate_fd_loop.m": loop + "\n% edited in the code view\n",
140+ });
141+ expect(g.compute?.loop).toBe("numbl");
142+ expect(g.compute?.note).toMatch(/edited/);
143+ expect(rg.info.extra.loop).toBe("numbl");
144+ expect(relL2(rc.probe.w, rg.probe.w)).toBeLessThan(1e-12);
145+ });
146+});
147+
148+describe("explicit finite differences without a GPU", () => {
149+ it("runs the whole bundle in numbl when the GPU is asked for but absent", async () => {
150+ const spec = testSpec();
151+ const r: RunResult = await cpu.run({
152+ bundle: buildRunBundle(spec, FD_METHOD, { gpu: true }),
153+ gpu: buildGpuStages(spec, FD_METHOD)!,
154+ });
155+ expect(r.ok).toBe(true);
156+ if (!r.ok) return;
157+ expect(r.compute?.loop).toBe("numbl");
158+ expect(r.compute?.note).toMatch(/not available/);
159+ const res = parsePlateResult(r.value);
160+ expect(res.info.extra.loop).toBe("numbl");
161+ expect(res.info.stable).toBe(true);
162+ });
163+
164+ it("writes the staged script and the plain script for the same problem", () => {
165+ const spec = testSpec();
166+ const plain = buildRunBundle(spec, FD_METHOD).files[0].source;
167+ const staged = buildRunBundle(spec, FD_METHOD, { gpu: true }).files[0].source;
168+ expect(plain).toMatch(/result = plate_fd\(params\);/);
169+ expect(staged).toMatch(/S = plate_fd_setup\(params\);/);
170+ expect(staged).toMatch(/R = plate_fd_loop\(S\);/);
171+ expect(staged).toMatch(/result = plate_fd_output\(S, R\);/);
172+ const stages = buildGpuStages(spec, FD_METHOD)!;
173+ expect(stages.loopFile).toBe("plate_fd_loop.m");
174+ expect(stages.setup).toMatch(/find\(S\.Km\)/);
175+ expect(stages.output).toMatch(/plate_fd_output\(S, R\)/);
176+ });
177+});
packages/core/test/implicit.test.tsmodified+24−5View file
@@ -7,7 +7,8 @@ describe("implicit finite differences", () => {
77 const check = `
88 for bc = {'simply-supported', 'clamped'}
99 nx = 19; ny = 13; dx = 0.3/(nx-1); dy = 0.2/(ny-1);
10- [L, B] = plate_matrices(nx, ny, dx, dy, bc{1});
10+ dom = plate_domain(struct('shape', 'rectangle', 'Lx', 0.3, 'Ly', 0.2), nx, ny);
11+ [L, B] = plate_matrices(dom, bc{1});
1112 w = zeros(ny, nx); w(2:end-1, 2:end-1) = reshape(sin(1:(ny-2)*(nx-2)), ny-2, nx-2);
1213 Lw = plate_laplacian(plate_pad(w, bc{1}), dx, dy);
1314 Bw = zeros(ny, nx);
@@ -19,7 +20,13 @@ describe("implicit finite differences", () => {
1920 end
2021 out = struct('err', err, 'sp', sp);
2122 `;
22- const run = await runner.run({ bundle: { files: [{ name: "c.m", source: check }, ...LIBRARY_FILES], entry: "c.m", resultVariable: "out" } });
23+ const run = await runner.run({
24+ bundle: {
25+ files: [{ name: "c.m", source: check }, ...LIBRARY_FILES],
26+ entry: "c.m",
27+ resultVariable: "out",
28+ },
29+ });
2330 expect(run.ok).toBe(true);
2431 if (!run.ok || !isPlainStruct(run.value)) return;
2532 const err = run.value.fields.err;
@@ -33,7 +40,13 @@ describe("implicit finite differences", () => {
3340
3441 it("agrees with the explicit scheme on the same grid and with the modal solution", async () => {
3542 const spec = testSpec();
36- spec.implicit = { nx: spec.fd.nx, ny: spec.fd.ny, duration: spec.fd.duration, steps_per_sample: 2, theta: 0.25 };
43+ spec.implicit = {
44+ nx: spec.fd.nx,
45+ ny: spec.fd.ny,
46+ duration: spec.fd.duration,
47+ steps_per_sample: 2,
48+ theta: 0.25,
49+ };
3750 spec.output.duration = spec.fd.duration;
3851 spec.output.grid = [spec.fd.nx, spec.fd.ny];
3952 const imp = await solve(spec, "implicit");
@@ -42,8 +55,14 @@ describe("implicit finite differences", () => {
4255 expect(imp.result.method).toBe("implicit");
4356 expect(imp.result.info.stable).toBe(true);
4457 expect(imp.result.info.nsteps).toBeLessThan(exp.result.info.nsteps);
45- const cmpE = await runner.run({ bundle: buildCompareBundle(), inputs: { A: exp.run.ok ? exp.run.value : 0, B: imp.run.ok ? imp.run.value : 0 } });
46- const cmpM = await runner.run({ bundle: buildCompareBundle(), inputs: { A: modal.run.ok ? modal.run.value : 0, B: imp.run.ok ? imp.run.value : 0 } });
58+ const cmpE = await runner.run({
59+ bundle: buildCompareBundle(),
60+ inputs: { A: exp.run.ok ? exp.run.value : 0, B: imp.run.ok ? imp.run.value : 0 },
61+ });
62+ const cmpM = await runner.run({
63+ bundle: buildCompareBundle(),
64+ inputs: { A: modal.run.ok ? modal.run.value : 0, B: imp.run.ok ? imp.run.value : 0 },
65+ });
4766 expect(cmpE.ok && cmpM.ok).toBe(true);
4867 if (!cmpE.ok || !cmpM.ok) return;
4968 const e = parseCompareResult(cmpE.value);
packages/core/test/wav.test.tsadded+39−0View file
@@ -0,0 +1,39 @@
1+import { describe, expect, it } from "vitest";
2+import { encodeWav } from "../src/audio/wav";
3+
4+const str = (b: Uint8Array, off: number, n: number) =>
5+ String.fromCharCode(...b.subarray(off, off + n));
6+
7+describe("encodeWav", () => {
8+ it("writes a well-formed 16-bit mono PCM header", () => {
9+ const b = encodeWav(new Float32Array(10), 44100);
10+ const v = new DataView(b.buffer);
11+ expect(b.length).toBe(44 + 20);
12+ expect(str(b, 0, 4)).toBe("RIFF");
13+ expect(v.getUint32(4, true)).toBe(b.length - 8);
14+ expect(str(b, 8, 4)).toBe("WAVE");
15+ expect(str(b, 12, 4)).toBe("fmt ");
16+ expect(v.getUint16(20, true)).toBe(1);
17+ expect(v.getUint16(22, true)).toBe(1);
18+ expect(v.getUint32(24, true)).toBe(44100);
19+ expect(v.getUint32(28, true)).toBe(88200);
20+ expect(v.getUint16(32, true)).toBe(2);
21+ expect(v.getUint16(34, true)).toBe(16);
22+ expect(str(b, 36, 4)).toBe("data");
23+ expect(v.getUint32(40, true)).toBe(20);
24+ });
25+
26+ it("quantizes, clips, and silences non-finite samples", () => {
27+ const b = encodeWav([0, 1, -1, 0.5, 2, -3, NaN, Infinity], 8000);
28+ const v = new DataView(b.buffer);
29+ const s = (i: number) => v.getInt16(44 + 2 * i, true);
30+ expect(s(0)).toBe(0);
31+ expect(s(1)).toBe(32767);
32+ expect(s(2)).toBe(-32768);
33+ expect(s(3)).toBe(16384);
34+ expect(s(4)).toBe(32767);
35+ expect(s(5)).toBe(-32768);
36+ expect(s(6)).toBe(0);
37+ expect(s(7)).toBe(0);
38+ });
39+});
tsconfig.base.jsonmodified+1−1View file
@@ -12,6 +12,6 @@
1212 "isolatedModules": true,
1313 "resolveJsonModule": true,
1414 "jsx": "react-jsx",
15- "types": ["vite/client"]
15+ "types": ["vite/client", "@webgpu/types"]
1616 }
1717 }