Clean rebuild of the steel sheet simulator
Von Karman plate with in-plane inertia, Marguerre bow, rectangular or
elliptical outline, push-and-release and drumstick excitations, and the
Rayleigh listener, stepped explicitly with WebGPU. Compared with
steelsim21: the drumstick contact uses the mass-weighted patch
displacement as its coordinate so that the contact has a potential
energy; the k^4 loss, the dipole heuristic, the command-line renderer,
the gallery, and the matching scripts are left out; the documentation is
rewritten with the grid dispersion stated explicitly; and a Playwright
test compares one GPU step with the float64 reference.
Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
31 changed files+5316−0
.github/workflows/deploy.ymladded+38−0View file
@@ -0,0 +1,38 @@
1+name: deploy
2+on:
3+ push:
4+ branches: [main]
5+ workflow_dispatch:
6+
7+permissions:
8+ contents: read
9+ pages: write
10+ id-token: write
11+
12+concurrency:
13+ group: pages
14+ cancel-in-progress: false
15+
16+jobs:
17+ build-deploy:
18+ runs-on: ubuntu-latest
19+ environment:
20+ name: github-pages
21+ url: ${{ steps.deployment.outputs.page_url }}
22+ steps:
23+ - uses: actions/checkout@v4
24+ - uses: actions/setup-node@v4
25+ with:
26+ node-version: 24
27+ cache: npm
28+ - run: npm ci
29+ # Discretization tests on the float64 reference (no GPU needed).
30+ - run: npm run test:cpu
31+ # Type check and Vite build.
32+ - run: npm run build
33+ - uses: actions/configure-pages@v5
34+ - uses: actions/upload-pages-artifact@v3
35+ with:
36+ path: dist
37+ - id: deployment
38+ uses: actions/deploy-pages@v4
.gitignoreadded+2−0View file
@@ -0,0 +1,2 @@
1+node_modules/
2+dist/
README.mdadded+49−0View file
@@ -0,0 +1,49 @@
1+# steelsim31
2+
3+A browser simulation of the sound of a thin stainless steel sheet, rectangular or elliptical, that hangs from a thread and is either pushed at a corner and released or struck with a drumstick. The sheet rings like a gong, so the simulator solves the geometrically nonlinear (von Karman) plate equations with in-plane inertia retained, stepped explicitly on the GPU through WebGPU. The model and the numerical method are described in [THEORY.md](THEORY.md).
4+
5+This is a clean rebuild of steelsim21, reduced to the parts that matter for the sound: the plate model with a stress-free bow, the two excitations, the two-parameter loss model plus an in-plane loss, the Rayleigh listener, and the tests. The sound-matching additions of steelsim21 (the k^4 loss term, the dipole heuristic, the matching scripts) and its command-line renderer and gallery were left out.
6+
7+The app is deployed by the `deploy` workflow on every push to `main`: https://steelsim31--proof-of-concept.vault1-sites.magland.org/
8+
9+## Running
10+
11+Requires Node 20 or later and a browser with WebGPU (Chrome or Edge 113+; Firefox and Safari with WebGPU enabled).
12+
13+ npm install
14+ npm run dev # development server
15+ npm run build # static site in dist/
16+ npm test # all tests (see below)
17+
18+## Using the simulator
19+
20+Press Start. With the push excitation the sheet is pushed at the bottom-right corner over the ramp time, held with heavy artificial damping until it settles, and released by ramping the finger force to zero over the release time; the recording of audio and movie frames starts at the release. With the drumstick excitation a free mass strikes the sheet at a configurable point and speed through a Hertz contact, and the recording starts at the impact. The view shows the transverse displacement as a color map (2D) or a lit height field (3D; drag to orbit, wheel to zoom), the spectrogram shows the listener signal, and the audio can be played or downloaded as a 16-bit WAV at any time.
21+
22+- The sheet is a rectangle by default; the Outline selector makes it the ellipse inscribed in the given width and height. Width and height are entered in inches by default, or in millimeters; everything else is metric. Internally all values are SI.
23+- The sheet has a stress-free bow (default 3 mm across the width; set the bow shape to Flat for the ideal plate). A negative push amount pushes against the bow.
24+- Sheet, initial shape, excitation, and numerics parameters take effect at the next Start. Losses and the listener position can be changed while a simulation runs.
25+- Pause halts the run temporarily; Stop ends it, keeps the audio and movie, and unlocks the parameter fields. To extend a finished run, raise the duration and press Continue; the state is kept on the GPU.
26+- Playback level: Absolute maps 1 Pa to full scale, shifted by the gain slider, so different runs can be compared by ear (the peak level in dBFS is shown, with a clipping warning); Loudness and Peak normalize each run instead. The WAV download uses the same setting.
27+- After the simulation completes, drag the time slider to explore the movie; Play starts the audio from the slider position and moves the movie and the spectrogram cursor with it. Live returns to the running view.
28+- Any parameter can be preset in the URL in SI units, for example `?spacing=0.006&duration=2&autostart=1` or `?excitation=strike&shape=ellipse`.
29+
30+The default grid spacing of 8 mm (39 by 58 nodes) keeps the simulation interactive but lowers the bending partials by 3 percent near 1 kHz and by 10 percent near 3.6 kHz through grid dispersion; the Sheet panel reports these two frequencies for the chosen spacing. The stable time step is set by the in-plane wave speed, about 1.3 microseconds at 8 mm, so 4 s of sound take about 3 million steps. On an Intel Iris Xe integrated GPU this runs at about 33k steps per second, or roughly 1.5 minutes for 4 s of sound; a discrete GPU is considerably faster. Halving the spacing to 4 mm costs about eight times as much (about 5 minutes for 4 s on the same GPU) and moves the 3 percent point to about 4 kHz.
31+
32+## Tests
33+
34+`npm test` runs two groups of tests with Node's built-in runner.
35+
36+- `tests/cpuRef.test.ts` checks the float64 reference of the spatial discretization without a GPU: symmetry of the stiffness with exactly three rigid modes per block, forces as the exact gradient of the discrete energy (flat and bowed, rectangle and ellipse), dissipativity of the loss terms, and convergence of the free square and circular plate frequencies to the reference values of Gorman and Leissa.
37+- `tests/gpu.test.ts` serves the app with Vite and drives a headless Chromium through Playwright. It compares one step of the WGSL kernel with the float64 reference on flat, bowed, and elliptical sheets (`?selftest=1` in the browser shows the same numbers), and runs a short push and a short strike to completion. It needs Playwright's Chromium (`npx playwright install chromium`) and a machine whose Chromium has WebGPU; `npm run test:cpu` runs only the first group.
38+
39+## Layout
40+
41+- `src/params.ts`: parameters, grid and mask, material constants, geometry of the excitation, grid dispersion estimate.
42+- `src/sim/cpuRef.ts`: float64 reference of the spatial discretization; used for the time-step estimate and the tests.
43+- `src/sim/shaders/`: WGSL kernels (time step, rigid-body projection, listener sum, display extraction).
44+- `src/sim/simulation.ts`: WebGPU driver (buffers, batch encoding, phases, readback).
45+- `src/selftest.ts`: the in-browser GPU-versus-CPU comparison.
46+- `src/render/`: 2D and 3D WebGPU rendering of the displacement field.
47+- `src/audio/`: FFT, spectrogram, WAV encoder, playback, level mapping.
48+- `src/ui/form.ts` and `src/main.ts`: parameter form and the rest of the interface.
49+- `tests/`: the two test files above.
THEORY.mdadded+87−0View file
@@ -0,0 +1,87 @@
1+# Model and numerical method
2+
3+This note describes the model behind the simulator: a thin stainless steel sheet, rectangular or elliptical, hanging from a thread, that is either pushed at a point of its rim, held, and released, or struck with a drumstick. The sheet rings like a gong. Its vibration amplitude, several thicknesses, is beyond the range where a linear plate model is accurate, so we solve the geometrically nonlinear plate equations of Foppl and von Karman, keep the in-plane inertia, discretize on a regular grid, and step explicitly on the GPU. Section 1 gives the continuous model, Sections 2 and 3 the suspension and the two excitations, Sections 4 and 5 the discretization in space and time, Section 6 the radiation model, and Section 7 the limitations.
4+
5+## 1. Plate model
6+
7+The sheet occupies the rectangle $[-L_x/2, L_x/2] \times [-L_y/2, L_y/2]$, or the ellipse inscribed in it, and has thickness $h$, Young's modulus $E$, Poisson ratio $\nu$, and density $\rho$. Its motion is described by the in-plane displacements $u(x,y,t)$ and $v(x,y,t)$ and the transverse displacement $w(x,y,t)$, all measured from a stress-free initial shape $w_0(x,y)$ that represents the bow of the real sheet (below). The membrane strains in the Marguerre form are
8+
9+$$
10+\varepsilon_{xx} = u_x + \tfrac12\left[(w_0 + w)_x^2 - w_{0,x}^2\right], \qquad
11+\varepsilon_{yy} = v_y + \tfrac12\left[(w_0 + w)_y^2 - w_{0,y}^2\right],
12+$$
13+
14+$$
15+\gamma_{xy} = u_y + v_x + (w_0 + w)_x (w_0 + w)_y - w_{0,x} w_{0,y},
16+$$
17+
18+so that $u = v = w = 0$ carries no stress. The membrane forces per unit length are $N_{xx} = C(\varepsilon_{xx} + \nu\varepsilon_{yy})$, $N_{yy} = C(\varepsilon_{yy} + \nu\varepsilon_{xx})$, and $N_{xy} = G\gamma_{xy}$, with $C = Eh/(1-\nu^2)$ and $G = Eh/(2(1+\nu))$. Writing $W = w_0 + w$ for the total shape, the equations of motion are
19+
20+$$
21+\rho h\, u_{tt} = \partial_x N_{xx} + \partial_y N_{xy}, \qquad
22+\rho h\, v_{tt} = \partial_x N_{xy} + \partial_y N_{yy},
23+$$
24+
25+$$
26+\rho h\, w_{tt} = -D\nabla^4 w + \partial_x(N_{xx} W_x + N_{xy} W_y) + \partial_y(N_{xy} W_x + N_{yy} W_y) - 2\sigma_0\rho h\, w_t + 2\sigma_1 \rho h\, \nabla^2 w_t ,
27+$$
28+
29+with the bending stiffness $D = Eh^3/(12(1-\nu^2))$. Bending acts on the deviation $w$ from the initial shape; the membrane forces act through the total slope. For $w_0 = 0$ these are the von Karman equations. The quadratic terms in the strains couple bending to stretching and produce the gong-like behavior: the pitch glide of a strongly excited plate and the transfer of energy from the low modes toward high frequencies. The in-plane inertia terms $\rho h u_{tt}$ and $\rho h v_{tt}$ are usually dropped in favor of an Airy stress function, since in-plane waves are much faster than bending waves. We keep them because they make the whole system explicit and local, which suits the GPU, at the cost of a time step set by the in-plane wave speed (Section 5).
30+
31+**Initial shape.** Real sheet metal is never flat; rolling and handling leave a bow of a few millimeters. A bowed sheet is a shallow shell, and its behavior differs qualitatively from a flat plate: the flat plate has a symmetric, purely hardening nonlinearity, whereas a shell has an asymmetric one, may snap between two curved states, and oscillates more slowly at large amplitude as it passes through the soft flat configuration. This is the mechanism behind oil-canning, and it appears to be what produces the slow wobble of a real sheet. Three parabolic shapes are offered: curved across the width, curved across the height, and a dish with the same curvature both ways. The sag is prescribed at the center of the bounding box; for the two cylindrical shapes it is measured relative to the straight edges, and for the dish relative to the corners, with the edge midpoints at half the sag. The default is a 3 mm bow across the width. A positive push moves the rim point in the direction of the bow and a negative push against it; the two cases sound different, as they do on a real sheet.
32+
33+**Losses.** The two-parameter loss model is the standard one in plate sound synthesis. A bending mode of angular frequency $\omega$ and wavenumber $k$ decays at the rate $\sigma_0 + \sigma_1 k^2$, and for the flat plate $k^2 = \omega/\kappa$ with $\kappa = \sqrt{D/\rho h}$, so $\sigma_0$ sets the low-frequency decay and $\sigma_1$ makes high frequencies decay faster in proportion to frequency. The defaults, $\sigma_0 = 1.15\ \mathrm{s^{-1}}$ and $\sigma_1 = 5\times10^{-4}\ \mathrm{m^2/s}$, give 60 dB decay times of about 4.9 s at 100 Hz, 1.8 s at 1 kHz, and 0.5 s at 5 kHz. The in-plane motion carries an additional frequency-independent loss $\sigma_{\mathrm{in}}$ (default $20\ \mathrm{s^{-1}}$). The stretching modes of the sheet lie at several kilohertz, couple linearly to the transverse motion once the sheet is curved, and would otherwise ring for the whole low-frequency decay time; a material loss factor of about $10^{-3}$ at 5 kHz gives a rate of this order. None of these values are measured; they are chosen to sound plausible, and all three are exposed in the interface.
34+
35+## 2. Suspension and rigid-body motion
36+
37+The whole edge is free: the bending moment and the Kirchhoff shear vanish, and so does the in-plane traction. Rather than imposing these conditions explicitly, the discretization of Section 4 is derived from a discrete energy, and free edges are then its natural boundary conditions.
38+
39+The thread through the hole near one corner is not modeled as a constraint. At audio frequencies it is a very soft spring compared with the local stiffness of the plate, so it does not pin the hole; it only fixes the average position of the plate. We therefore treat the plate as free and remove the six rigid-body modes (three translations, three rotations) from the displacement and velocity fields by a mass-weighted projection every few time steps. Gravity is dropped: the in-plane stress it induces is about $\rho g L \approx 35$ kPa, a strain of $2\times10^{-7}$.
40+
41+The same projection defines what holding the sheet means during a push. The finger prescribes the displacement of a small patch, and the reaction to the finger force is distributed over the plate as an inertial load, which is the inertia-relief equilibrium of a free body, so the push produces pure bending of the prescribed amount. This is a modeling choice: a sheet hanging freely from a thread and pushed 10 mm at a corner would mostly swing as a pendulum and bend by only a millimeter or two. The same choice applies to the drumstick, whose force would otherwise set the whole sheet swinging.
42+
43+## 3. Excitation
44+
45+**Push and release.** The simulation has three phases. During the push (default 0.1 s) the transverse displacement of the nodes within the finger radius follows a smooth ramp to the push amount, and heavy artificial damping ($\sigma_0 = 300\ \mathrm{s^{-1}}$) removes the waves that the ramp launches. During the hold the constraint stays in place until the non-rigid kinetic energy drops below a threshold or a maximum settle time elapses. At the release the clock is set to zero, the damping returns to its physical value, and the recording of audio and movie frames begins. The finger does not let go instantaneously: the reaction force that held each node of the patch is recorded at the end of the hold and then applied as an external force that falls smoothly to zero over the release time (20 ms by default). This matters for the sound. An instantaneous release is a step in the finger force, which excites all modes at once, and since the radiated pressure is proportional to the plate acceleration, the result is a nearly flat pressure spectrum, a cymbal-like click. A release over 20 ms confines the excitation to roughly the first hundred hertz, as a fingertip does. The rigid-body components of the recorded reaction are removed before it is applied, so that the finger force and the distributed holding reaction relax together. The release time can be set to zero to recover the instantaneous case.
46+
47+**Drumstick strike.** The stick tip is a free mass $M$ (default 30 g, the effective mass of a stick in a wrist stroke rather than its full mass) arriving with speed $v_0$ (default 3 m/s) at a point of the sheet. The contact patch consists of the sheet nodes within the tip radius (5 mm) plus the nearest node. Its contact coordinate is the mass-weighted mean transverse displacement $w_c$ of these nodes, the penetration is $\delta = w_c - x$ with $x$ the tip position, and the force is the Hertz law $F = K\,\delta^{\alpha}$ for $\delta > 0$ and zero otherwise, with $\alpha = 1.5$ for a rounded tip and $K = 5\times10^{8}\ \mathrm{N/m^{1.5}}$ for wood on steel. The force is spread over the patch in proportion to the node masses. Because that distribution is exactly conjugate to the mean displacement, the contact has a potential energy and the coupled system is Hamiltonian apart from the losses. The stick obeys $M\ddot x = F$ and is stepped together with the plate; it may bounce and strike again if the sheet comes back to meet it. The recording starts at the impact, with the sheet at rest in its natural shape. Note that a thin sheet is a soft target: its driving-point impedance is $8\sqrt{D\rho h} \approx 59\ \mathrm{N\,s/m}$, so the contact lasts well under a millisecond and the peak force is of the order of a hundred newtons. The excitation is therefore broadband, unlike the push, and the sound is a short bright attack followed by the ring.
48+
49+## 4. Spatial discretization
50+
51+The plate is covered by a regular grid of $n_x \times n_y$ nodes with spacings $h_x$ and $h_y$ (8 mm by default). All three displacements live at the nodes. The forces are the negative gradient of a discrete potential energy $U = U_b + U_m$, so the semi-discrete system is Hamiltonian, its linear part has a symmetric stiffness matrix, and the free edge needs no special formulas.
52+
53+**Bending.** The continuous bending energy is $\tfrac{D}{2}\int [w_{xx}^2 + w_{yy}^2 + 2\nu w_{xx}w_{yy} + 2(1-\nu) w_{xy}^2]\,dA$. We evaluate $w_{xx}$ and $w_{yy}$ by three-point second differences at the nodes where both neighbors exist, and $w_{xy}$ by the four-point difference at cell centers, with trapezoidal quadrature weights (1 inside, 1/2 on an edge, 1/4 at a corner). Where a second difference does not exist, at the edges perpendicular to the direction of differencing, its term is simply absent from the energy. In the interior the resulting force is the standard 13-point biharmonic stencil; near the edges the stencil is modified exactly as the energy dictates. The three rigid modes $w = 1, x, y$ lie in the null space.
54+
55+**Membrane.** The cells are bilinear (Q4) finite elements with 2 by 2 Gauss integration, using the strains of Section 1 with $w_0$ and $w$ interpolated bilinearly. One-point integration would be cheaper but admits hourglass modes, which the nonlinear coupling would excite. On the GPU the force of each cell on its four nodes is computed once per cell within a workgroup and gathered by the nodes.
56+
57+**Mass.** The mass matrix is lumped with the same trapezoidal weights, so that the total mass is exact.
58+
59+**Elliptical sheets.** The same grid serves an elliptical sheet through a mask. A cell belongs to the sheet when its center lies inside the ellipse, and a node when it is a corner of such a cell; the edge is therefore a staircase of whole cells. Every rule above carries over unchanged: a second difference exists where the node and both of its neighbors belong to the sheet, a cell contributes membrane energy and the mixed bending term when its four corners do, and the weight of a node is the fraction of its four cells that belong to the sheet. Because the forces are still the gradient of the discrete energy, the staircase edge is free, the stiffness stays symmetric, and exactly three rigid modes remain in each block. Nodes outside the sheet are frozen and carry no mass. Deciding membership by cell centers rather than by nodes keeps the area unbiased; a node-based mask loses about half a cell all around the edge and raises every frequency by several percent. The hole and the pushed patch are placed on the rim in the upper-left and lower-right directions, so the push remains asymmetric as it is at the corner of a rectangle.
60+
61+**Accuracy.** Two errors limit the accuracy of the spatial scheme, and it is worth separating them. The first is the treatment of the edge. For the rectangle the scheme converges at second order: the first seven frequencies of a completely free square plate are low by about 3.5 percent at 10 cells per side, 0.8 percent at 20, and 0.3 percent at 30, compared with the reference values of Gorman. The staircase edge of the ellipse is only first-order accurate; the frequencies of a free circular plate are within about 2 percent of Leissa's values at 24 cells across the diameter and within 1 percent at 60, and the degenerate pairs, which the staircase splits, differ by 1 to 2 percent. The second error is grid dispersion, which affects every partial regardless of the outline. A plane bending wave of wavenumber $k$ along a grid axis has the discrete frequency ratio $(2 - 2\cos kh)/(kh)^2$, so a partial comes out low by 1 percent at $kh \approx 0.35$, by 3 percent at $kh \approx 0.60$, and by 10 percent at $kh \approx 1.12$. For the default sheet this places the 3 percent point near 1.1 kHz on the 8 mm grid, 1.9 kHz on a 6 mm grid, and 4.2 kHz on a 4 mm grid, and the 10 percent point near 3.6, 6.5, and 15 kHz respectively. The high partials of a gong sound are therefore audibly flattened at 8 mm, and the interface reports these two frequencies for the chosen spacing. Since the cost scales with the inverse cube of the spacing (Section 5), the default is a compromise in favor of interactive use, and a finer grid is a matter of patience rather than of the method.
62+
63+**Validation.** The tests check that both stiffness blocks are symmetric with exactly three rigid modes, that the forces are the exact gradient of the discrete energy at a strongly nonlinear state with and without a bow, that the loss terms only remove energy, that the free-plate frequencies converge as stated, and that the same holds on an elliptical mask. A further test runs one time step of the GPU kernel on a smooth random state of about 2.5 thicknesses amplitude, flat, bowed, and elliptical, and compares it with the float64 reference: the accelerations agree to about $10^{-6}$ relative, which is single-precision rounding.
64+
65+## 5. Time integration
66+
67+The system is advanced with the symplectic Euler scheme, $v^{n+1} = v^n + \Delta t\, a(x^n, v^n)$ and $x^{n+1} = x^n + \Delta t\, v^{n+1}$, in single precision on the GPU. With the velocities interpreted at half steps this is the leapfrog (Stormer-Verlet) scheme, second order and symplectic for the conservative part; the losses are evaluated explicitly at the old velocity. The stability limit is $\Delta t < 2/\omega_{\max}$, where $\omega_{\max}$ is the largest frequency of the linearized discrete system. We estimate $\omega_{\max}$ at start-up by power iteration on a small grid with the same spacing and use 85 percent of the limit. The in-plane waves, with speed $c_p = \sqrt{E/\rho(1-\nu^2)} \approx 5200$ m/s, set this limit: at 8 mm spacing $\Delta t \approx 1.3\ \mu$s, roughly six times smaller than bending alone would require, and 4 s of sound take about 3 million steps; at 4 mm the step halves and the node count quadruples. The step is rounded down so that an integer number of steps (16 at 8 mm and 48 kHz) makes one audio sample.
68+
69+Single precision needs one precaution. With $w$ of a few millimeters and $\Delta t \approx 1\ \mu$s, a velocity below about $10^{-3}$ m/s produces an increment smaller than the resolution of $w$, and the rounding would act as a broadband noise source. The update of $w$ is therefore compensated: the rounding error of each addition is carried in a separate array and added back at the next step (Kahan summation). The self-test of Section 4 also checks this update.
70+
71+The state is stored twice and each step reads one copy and writes the other. A step is one compute dispatch of 16 by 16 workgroups; a per-step uniform, addressed by a dynamic offset, carries the prescribed patch displacement, the damping, and the audio bookkeeping, so that a whole movie frame of steps is encoded in one pass.
72+
73+## 6. Radiation to the listener
74+
75+The listener model is the Rayleigh integral for a baffled plate,
76+
77+$$
78+p(\mathbf r, t) = \frac{\rho_{\mathrm{air}}}{2\pi} \sum_k \frac{\Delta A_k}{r_k}\, \ddot w_k(t - r_k/c),
79+$$
80+
81+where the sum runs over the nodes, $r_k$ is the distance from node $k$ to the listener, and $c = 343$ m/s. The transverse velocity of every node is recorded in a ring buffer at the audio rate, and each output sample is a delayed, weighted sum over the plate, evaluated on the GPU with linear interpolation between samples; the acceleration is formed as the difference between successive delayed velocities, which is the average acceleration over one sample interval. This matters in single precision: the force evaluation carries a broadband rounding error of the order of $1\ \mathrm{m/s^2}$ per node, and sampling the instantaneous acceleration would put that noise straight into the audio, whereas the averaged form suppresses it by the square root of the number of steps per sample. The common delay is discarded; only the spread of delays across the plate (about 5 samples for a listener 1 m in front) is kept.
82+
83+The plate is unbaffled in reality, so the sound from the two faces partly cancels at low frequencies, and the Rayleigh integral overestimates the low-frequency level. An accurate treatment would require solving the exterior acoustic problem, which is outside the scope of this simulator. The pressure is reported in pascals and the peak sound pressure level is displayed. Playback and the WAV file use, by default, a fixed mapping from pressure to full scale (1 Pa at 0 dBFS, shifted by a user gain in decibels), so that runs with different excitations can be compared by ear; a soft clipper bounds the output and the resulting peak level is shown. Two normalized modes are also offered, to the peak or to the loudness of the sustained sound.
84+
85+## 7. Limitations
86+
87+The von Karman model assumes moderate rotations, so it is not valid for displacements that are large compared with the plate size; a corner push of 10 mm on a 300 by 460 mm sheet is well within its range. The Marguerre treatment of the initial shape is a shallow-shell approximation, adequate for a bow of a few millimeters, and the real shape of a given sheet is of course unknown; the three offered shapes are idealizations. The loss model is empirical and its parameters are not calibrated against a real sheet. Grid dispersion lowers the high partials by the amounts given in Section 4, and the staircase edge of an elliptical sheet places the low frequencies within a few percent and splits modes that would be degenerate on a true ellipse. The suspension is idealized as inertia relief, so the swinging of the sheet on its thread is absent. Air loading, which lowers the frequencies of a thin plate slightly and adds radiation damping, is not included, and the radiation model is approximate, as discussed above. Finally, single precision limits the dynamic range of the tail of the sound to roughly 100 dB below the peak.
index.htmladded+29−0View file
@@ -0,0 +1,29 @@
1+<!doctype html>
2+<html lang="en">
3+ <head>
4+ <meta charset="UTF-8" />
5+ <meta name="viewport" content="width=device-width, initial-scale=1.0" />
6+ <title>Steel sheet simulator</title>
7+ <link rel="stylesheet" href="./src/style.css" />
8+ </head>
9+ <body>
10+ <div id="app">
11+ <aside id="sidebar">
12+ <h1>Steel sheet</h1>
13+ <p class="sub">Geometrically nonlinear plate (von Karman with in-plane inertia), stepped explicitly on the GPU.</p>
14+ <div id="controls"></div>
15+ <div id="params"></div>
16+ <div id="info"></div>
17+ <p class="sub">The model and the numerical method are described in THEORY.md in the repository.</p>
18+ </aside>
19+ <main id="main">
20+ <div class="toolbar" id="viewToolbar"></div>
21+ <div id="viewWrap"><canvas id="view"></canvas><div id="nogpu" hidden></div></div>
22+ <div class="toolbar" id="movieToolbar"></div>
23+ <div class="toolbar" id="specToolbar"></div>
24+ <canvas id="spec"></canvas>
25+ </main>
26+ </div>
27+ <script type="module" src="./src/main.ts"></script>
28+ </body>
29+</html>
package-lock.jsonadded+1235−0View file
@@ -0,0 +1,1235 @@
1+{
2+ "name": "steelsim31",
3+ "version": "0.1.0",
4+ "lockfileVersion": 3,
5+ "requires": true,
6+ "packages": {
7+ "": {
8+ "name": "steelsim31",
9+ "version": "0.1.0",
10+ "devDependencies": {
11+ "@types/node": "^24.0.0",
12+ "@webgpu/types": "^0.1.60",
13+ "playwright": "^1.50.0",
14+ "typescript": "^5.6.0",
15+ "vite": "^6.0.0"
16+ }
17+ },
18+ "node_modules/@esbuild/aix-ppc64": {
19+ "version": "0.25.12",
20+ "resolved": "https://registry.npmjs.org/@esbuild/aix-ppc64/-/aix-ppc64-0.25.12.tgz",
21+ "integrity": "sha512-Hhmwd6CInZ3dwpuGTF8fJG6yoWmsToE+vYgD4nytZVxcu1ulHpUQRAB1UJ8+N1Am3Mz4+xOByoQoSZf4D+CpkA==",
22+ "cpu": [
23+ "ppc64"
24+ ],
25+ "dev": true,
26+ "license": "MIT",
27+ "optional": true,
28+ "os": [
29+ "aix"
30+ ],
31+ "engines": {
32+ "node": ">=18"
33+ }
34+ },
35+ "node_modules/@esbuild/android-arm": {
36+ "version": "0.25.12",
37+ "resolved": "https://registry.npmjs.org/@esbuild/android-arm/-/android-arm-0.25.12.tgz",
38+ "integrity": "sha512-VJ+sKvNA/GE7Ccacc9Cha7bpS8nyzVv0jdVgwNDaR4gDMC/2TTRc33Ip8qrNYUcpkOHUT5OZ0bUcNNVZQ9RLlg==",
39+ "cpu": [
40+ "arm"
41+ ],
42+ "dev": true,
43+ "license": "MIT",
44+ "optional": true,
45+ "os": [
46+ "android"
47+ ],
48+ "engines": {
49+ "node": ">=18"
50+ }
51+ },
52+ "node_modules/@esbuild/android-arm64": {
53+ "version": "0.25.12",
54+ "resolved": "https://registry.npmjs.org/@esbuild/android-arm64/-/android-arm64-0.25.12.tgz",
55+ "integrity": "sha512-6AAmLG7zwD1Z159jCKPvAxZd4y/VTO0VkprYy+3N2FtJ8+BQWFXU+OxARIwA46c5tdD9SsKGZ/1ocqBS/gAKHg==",
56+ "cpu": [
57+ "arm64"
58+ ],
59+ "dev": true,
60+ "license": "MIT",
61+ "optional": true,
62+ "os": [
63+ "android"
64+ ],
65+ "engines": {
66+ "node": ">=18"
67+ }
68+ },
69+ "node_modules/@esbuild/android-x64": {
70+ "version": "0.25.12",
71+ "resolved": "https://registry.npmjs.org/@esbuild/android-x64/-/android-x64-0.25.12.tgz",
72+ "integrity": "sha512-5jbb+2hhDHx5phYR2By8GTWEzn6I9UqR11Kwf22iKbNpYrsmRB18aX/9ivc5cabcUiAT/wM+YIZ6SG9QO6a8kg==",
73+ "cpu": [
74+ "x64"
75+ ],
76+ "dev": true,
77+ "license": "MIT",
78+ "optional": true,
79+ "os": [
80+ "android"
81+ ],
82+ "engines": {
83+ "node": ">=18"
84+ }
85+ },
86+ "node_modules/@esbuild/darwin-arm64": {
87+ "version": "0.25.12",
88+ "resolved": "https://registry.npmjs.org/@esbuild/darwin-arm64/-/darwin-arm64-0.25.12.tgz",
89+ "integrity": "sha512-N3zl+lxHCifgIlcMUP5016ESkeQjLj/959RxxNYIthIg+CQHInujFuXeWbWMgnTo4cp5XVHqFPmpyu9J65C1Yg==",
90+ "cpu": [
91+ "arm64"
92+ ],
93+ "dev": true,
94+ "license": "MIT",
95+ "optional": true,
96+ "os": [
97+ "darwin"
98+ ],
99+ "engines": {
100+ "node": ">=18"
101+ }
102+ },
103+ "node_modules/@esbuild/darwin-x64": {
104+ "version": "0.25.12",
105+ "resolved": "https://registry.npmjs.org/@esbuild/darwin-x64/-/darwin-x64-0.25.12.tgz",
106+ "integrity": "sha512-HQ9ka4Kx21qHXwtlTUVbKJOAnmG1ipXhdWTmNXiPzPfWKpXqASVcWdnf2bnL73wgjNrFXAa3yYvBSd9pzfEIpA==",
107+ "cpu": [
108+ "x64"
109+ ],
110+ "dev": true,
111+ "license": "MIT",
112+ "optional": true,
113+ "os": [
114+ "darwin"
115+ ],
116+ "engines": {
117+ "node": ">=18"
118+ }
119+ },
120+ "node_modules/@esbuild/freebsd-arm64": {
121+ "version": "0.25.12",
122+ "resolved": "https://registry.npmjs.org/@esbuild/freebsd-arm64/-/freebsd-arm64-0.25.12.tgz",
123+ "integrity": "sha512-gA0Bx759+7Jve03K1S0vkOu5Lg/85dou3EseOGUes8flVOGxbhDDh/iZaoek11Y8mtyKPGF3vP8XhnkDEAmzeg==",
124+ "cpu": [
125+ "arm64"
126+ ],
127+ "dev": true,
128+ "license": "MIT",
129+ "optional": true,
130+ "os": [
131+ "freebsd"
132+ ],
133+ "engines": {
134+ "node": ">=18"
135+ }
136+ },
137+ "node_modules/@esbuild/freebsd-x64": {
138+ "version": "0.25.12",
139+ "resolved": "https://registry.npmjs.org/@esbuild/freebsd-x64/-/freebsd-x64-0.25.12.tgz",
140+ "integrity": "sha512-TGbO26Yw2xsHzxtbVFGEXBFH0FRAP7gtcPE7P5yP7wGy7cXK2oO7RyOhL5NLiqTlBh47XhmIUXuGciXEqYFfBQ==",
141+ "cpu": [
142+ "x64"
143+ ],
144+ "dev": true,
145+ "license": "MIT",
146+ "optional": true,
147+ "os": [
148+ "freebsd"
149+ ],
150+ "engines": {
151+ "node": ">=18"
152+ }
153+ },
154+ "node_modules/@esbuild/linux-arm": {
155+ "version": "0.25.12",
156+ "resolved": "https://registry.npmjs.org/@esbuild/linux-arm/-/linux-arm-0.25.12.tgz",
157+ "integrity": "sha512-lPDGyC1JPDou8kGcywY0YILzWlhhnRjdof3UlcoqYmS9El818LLfJJc3PXXgZHrHCAKs/Z2SeZtDJr5MrkxtOw==",
158+ "cpu": [
159+ "arm"
160+ ],
161+ "dev": true,
162+ "license": "MIT",
163+ "optional": true,
164+ "os": [
165+ "linux"
166+ ],
167+ "engines": {
168+ "node": ">=18"
169+ }
170+ },
171+ "node_modules/@esbuild/linux-arm64": {
172+ "version": "0.25.12",
173+ "resolved": "https://registry.npmjs.org/@esbuild/linux-arm64/-/linux-arm64-0.25.12.tgz",
174+ "integrity": "sha512-8bwX7a8FghIgrupcxb4aUmYDLp8pX06rGh5HqDT7bB+8Rdells6mHvrFHHW2JAOPZUbnjUpKTLg6ECyzvas2AQ==",
175+ "cpu": [
176+ "arm64"
177+ ],
178+ "dev": true,
179+ "license": "MIT",
180+ "optional": true,
181+ "os": [
182+ "linux"
183+ ],
184+ "engines": {
185+ "node": ">=18"
186+ }
187+ },
188+ "node_modules/@esbuild/linux-ia32": {
189+ "version": "0.25.12",
190+ "resolved": "https://registry.npmjs.org/@esbuild/linux-ia32/-/linux-ia32-0.25.12.tgz",
191+ "integrity": "sha512-0y9KrdVnbMM2/vG8KfU0byhUN+EFCny9+8g202gYqSSVMonbsCfLjUO+rCci7pM0WBEtz+oK/PIwHkzxkyharA==",
192+ "cpu": [
193+ "ia32"
194+ ],
195+ "dev": true,
196+ "license": "MIT",
197+ "optional": true,
198+ "os": [
199+ "linux"
200+ ],
201+ "engines": {
202+ "node": ">=18"
203+ }
204+ },
205+ "node_modules/@esbuild/linux-loong64": {
206+ "version": "0.25.12",
207+ "resolved": "https://registry.npmjs.org/@esbuild/linux-loong64/-/linux-loong64-0.25.12.tgz",
208+ "integrity": "sha512-h///Lr5a9rib/v1GGqXVGzjL4TMvVTv+s1DPoxQdz7l/AYv6LDSxdIwzxkrPW438oUXiDtwM10o9PmwS/6Z0Ng==",
209+ "cpu": [
210+ "loong64"
211+ ],
212+ "dev": true,
213+ "license": "MIT",
214+ "optional": true,
215+ "os": [
216+ "linux"
217+ ],
218+ "engines": {
219+ "node": ">=18"
220+ }
221+ },
222+ "node_modules/@esbuild/linux-mips64el": {
223+ "version": "0.25.12",
224+ "resolved": "https://registry.npmjs.org/@esbuild/linux-mips64el/-/linux-mips64el-0.25.12.tgz",
225+ "integrity": "sha512-iyRrM1Pzy9GFMDLsXn1iHUm18nhKnNMWscjmp4+hpafcZjrr2WbT//d20xaGljXDBYHqRcl8HnxbX6uaA/eGVw==",
226+ "cpu": [
227+ "mips64el"
228+ ],
229+ "dev": true,
230+ "license": "MIT",
231+ "optional": true,
232+ "os": [
233+ "linux"
234+ ],
235+ "engines": {
236+ "node": ">=18"
237+ }
238+ },
239+ "node_modules/@esbuild/linux-ppc64": {
240+ "version": "0.25.12",
241+ "resolved": "https://registry.npmjs.org/@esbuild/linux-ppc64/-/linux-ppc64-0.25.12.tgz",
242+ "integrity": "sha512-9meM/lRXxMi5PSUqEXRCtVjEZBGwB7P/D4yT8UG/mwIdze2aV4Vo6U5gD3+RsoHXKkHCfSxZKzmDssVlRj1QQA==",
243+ "cpu": [
244+ "ppc64"
245+ ],
246+ "dev": true,
247+ "license": "MIT",
248+ "optional": true,
249+ "os": [
250+ "linux"
251+ ],
252+ "engines": {
253+ "node": ">=18"
254+ }
255+ },
256+ "node_modules/@esbuild/linux-riscv64": {
257+ "version": "0.25.12",
258+ "resolved": "https://registry.npmjs.org/@esbuild/linux-riscv64/-/linux-riscv64-0.25.12.tgz",
259+ "integrity": "sha512-Zr7KR4hgKUpWAwb1f3o5ygT04MzqVrGEGXGLnj15YQDJErYu/BGg+wmFlIDOdJp0PmB0lLvxFIOXZgFRrdjR0w==",
260+ "cpu": [
261+ "riscv64"
262+ ],
263+ "dev": true,
264+ "license": "MIT",
265+ "optional": true,
266+ "os": [
267+ "linux"
268+ ],
269+ "engines": {
270+ "node": ">=18"
271+ }
272+ },
273+ "node_modules/@esbuild/linux-s390x": {
274+ "version": "0.25.12",
275+ "resolved": "https://registry.npmjs.org/@esbuild/linux-s390x/-/linux-s390x-0.25.12.tgz",
276+ "integrity": "sha512-MsKncOcgTNvdtiISc/jZs/Zf8d0cl/t3gYWX8J9ubBnVOwlk65UIEEvgBORTiljloIWnBzLs4qhzPkJcitIzIg==",
277+ "cpu": [
278+ "s390x"
279+ ],
280+ "dev": true,
281+ "license": "MIT",
282+ "optional": true,
283+ "os": [
284+ "linux"
285+ ],
286+ "engines": {
287+ "node": ">=18"
288+ }
289+ },
290+ "node_modules/@esbuild/linux-x64": {
291+ "version": "0.25.12",
292+ "resolved": "https://registry.npmjs.org/@esbuild/linux-x64/-/linux-x64-0.25.12.tgz",
293+ "integrity": "sha512-uqZMTLr/zR/ed4jIGnwSLkaHmPjOjJvnm6TVVitAa08SLS9Z0VM8wIRx7gWbJB5/J54YuIMInDquWyYvQLZkgw==",
294+ "cpu": [
295+ "x64"
296+ ],
297+ "dev": true,
298+ "license": "MIT",
299+ "optional": true,
300+ "os": [
301+ "linux"
302+ ],
303+ "engines": {
304+ "node": ">=18"
305+ }
306+ },
307+ "node_modules/@esbuild/netbsd-arm64": {
308+ "version": "0.25.12",
309+ "resolved": "https://registry.npmjs.org/@esbuild/netbsd-arm64/-/netbsd-arm64-0.25.12.tgz",
310+ "integrity": "sha512-xXwcTq4GhRM7J9A8Gv5boanHhRa/Q9KLVmcyXHCTaM4wKfIpWkdXiMog/KsnxzJ0A1+nD+zoecuzqPmCRyBGjg==",
311+ "cpu": [
312+ "arm64"
313+ ],
314+ "dev": true,
315+ "license": "MIT",
316+ "optional": true,
317+ "os": [
318+ "netbsd"
319+ ],
320+ "engines": {
321+ "node": ">=18"
322+ }
323+ },
324+ "node_modules/@esbuild/netbsd-x64": {
325+ "version": "0.25.12",
326+ "resolved": "https://registry.npmjs.org/@esbuild/netbsd-x64/-/netbsd-x64-0.25.12.tgz",
327+ "integrity": "sha512-Ld5pTlzPy3YwGec4OuHh1aCVCRvOXdH8DgRjfDy/oumVovmuSzWfnSJg+VtakB9Cm0gxNO9BzWkj6mtO1FMXkQ==",
328+ "cpu": [
329+ "x64"
330+ ],
331+ "dev": true,
332+ "license": "MIT",
333+ "optional": true,
334+ "os": [
335+ "netbsd"
336+ ],
337+ "engines": {
338+ "node": ">=18"
339+ }
340+ },
341+ "node_modules/@esbuild/openbsd-arm64": {
342+ "version": "0.25.12",
343+ "resolved": "https://registry.npmjs.org/@esbuild/openbsd-arm64/-/openbsd-arm64-0.25.12.tgz",
344+ "integrity": "sha512-fF96T6KsBo/pkQI950FARU9apGNTSlZGsv1jZBAlcLL1MLjLNIWPBkj5NlSz8aAzYKg+eNqknrUJ24QBybeR5A==",
345+ "cpu": [
346+ "arm64"
347+ ],
348+ "dev": true,
349+ "license": "MIT",
350+ "optional": true,
351+ "os": [
352+ "openbsd"
353+ ],
354+ "engines": {
355+ "node": ">=18"
356+ }
357+ },
358+ "node_modules/@esbuild/openbsd-x64": {
359+ "version": "0.25.12",
360+ "resolved": "https://registry.npmjs.org/@esbuild/openbsd-x64/-/openbsd-x64-0.25.12.tgz",
361+ "integrity": "sha512-MZyXUkZHjQxUvzK7rN8DJ3SRmrVrke8ZyRusHlP+kuwqTcfWLyqMOE3sScPPyeIXN/mDJIfGXvcMqCgYKekoQw==",
362+ "cpu": [
363+ "x64"
364+ ],
365+ "dev": true,
366+ "license": "MIT",
367+ "optional": true,
368+ "os": [
369+ "openbsd"
370+ ],
371+ "engines": {
372+ "node": ">=18"
373+ }
374+ },
375+ "node_modules/@esbuild/openharmony-arm64": {
376+ "version": "0.25.12",
377+ "resolved": "https://registry.npmjs.org/@esbuild/openharmony-arm64/-/openharmony-arm64-0.25.12.tgz",
378+ "integrity": "sha512-rm0YWsqUSRrjncSXGA7Zv78Nbnw4XL6/dzr20cyrQf7ZmRcsovpcRBdhD43Nuk3y7XIoW2OxMVvwuRvk9XdASg==",
379+ "cpu": [
380+ "arm64"
381+ ],
382+ "dev": true,
383+ "license": "MIT",
384+ "optional": true,
385+ "os": [
386+ "openharmony"
387+ ],
388+ "engines": {
389+ "node": ">=18"
390+ }
391+ },
392+ "node_modules/@esbuild/sunos-x64": {
393+ "version": "0.25.12",
394+ "resolved": "https://registry.npmjs.org/@esbuild/sunos-x64/-/sunos-x64-0.25.12.tgz",
395+ "integrity": "sha512-3wGSCDyuTHQUzt0nV7bocDy72r2lI33QL3gkDNGkod22EsYl04sMf0qLb8luNKTOmgF/eDEDP5BFNwoBKH441w==",
396+ "cpu": [
397+ "x64"
398+ ],
399+ "dev": true,
400+ "license": "MIT",
401+ "optional": true,
402+ "os": [
403+ "sunos"
404+ ],
405+ "engines": {
406+ "node": ">=18"
407+ }
408+ },
409+ "node_modules/@esbuild/win32-arm64": {
410+ "version": "0.25.12",
411+ "resolved": "https://registry.npmjs.org/@esbuild/win32-arm64/-/win32-arm64-0.25.12.tgz",
412+ "integrity": "sha512-rMmLrur64A7+DKlnSuwqUdRKyd3UE7oPJZmnljqEptesKM8wx9J8gx5u0+9Pq0fQQW8vqeKebwNXdfOyP+8Bsg==",
413+ "cpu": [
414+ "arm64"
415+ ],
416+ "dev": true,
417+ "license": "MIT",
418+ "optional": true,
419+ "os": [
420+ "win32"
421+ ],
422+ "engines": {
423+ "node": ">=18"
424+ }
425+ },
426+ "node_modules/@esbuild/win32-ia32": {
427+ "version": "0.25.12",
428+ "resolved": "https://registry.npmjs.org/@esbuild/win32-ia32/-/win32-ia32-0.25.12.tgz",
429+ "integrity": "sha512-HkqnmmBoCbCwxUKKNPBixiWDGCpQGVsrQfJoVGYLPT41XWF8lHuE5N6WhVia2n4o5QK5M4tYr21827fNhi4byQ==",
430+ "cpu": [
431+ "ia32"
432+ ],
433+ "dev": true,
434+ "license": "MIT",
435+ "optional": true,
436+ "os": [
437+ "win32"
438+ ],
439+ "engines": {
440+ "node": ">=18"
441+ }
442+ },
443+ "node_modules/@esbuild/win32-x64": {
444+ "version": "0.25.12",
445+ "resolved": "https://registry.npmjs.org/@esbuild/win32-x64/-/win32-x64-0.25.12.tgz",
446+ "integrity": "sha512-alJC0uCZpTFrSL0CCDjcgleBXPnCrEAhTBILpeAp7M/OFgoqtAetfBzX0xM00MUsVVPpVjlPuMbREqnZCXaTnA==",
447+ "cpu": [
448+ "x64"
449+ ],
450+ "dev": true,
451+ "license": "MIT",
452+ "optional": true,
453+ "os": [
454+ "win32"
455+ ],
456+ "engines": {
457+ "node": ">=18"
458+ }
459+ },
460+ "node_modules/@napi-rs/lzma-linux-x64-gnu": {
461+ "version": "1.5.1",
462+ "resolved": "https://registry.npmjs.org/@napi-rs/lzma-linux-x64-gnu/-/lzma-linux-x64-gnu-1.5.1.tgz",
463+ "integrity": "sha512-oTXEIha4SsuXdTA4Iyskj0kpdx2yVXdhd75c2v3xGrHFfVMsbhTPZU/nMPL4sWKo4pBHm3aucLaqGlF696dTyQ==",
464+ "cpu": [
465+ "x64"
466+ ],
467+ "dev": true,
468+ "libc": [
469+ "glibc"
470+ ],
471+ "license": "MIT",
472+ "optional": true,
473+ "os": [
474+ "linux"
475+ ],
476+ "engines": {
477+ "node": "^22.20 || ^24.12 || >=25"
478+ }
479+ },
480+ "node_modules/@rollup/rollup-android-arm-eabi": {
481+ "version": "4.63.1",
482+ "resolved": "https://registry.npmjs.org/@rollup/rollup-android-arm-eabi/-/rollup-android-arm-eabi-4.63.1.tgz",
483+ "integrity": "sha512-UZ8sUxPTiHWYX9QNdJedb1kDZSpS1t/VPWBWGSgqHNi9w3Cu6IXvu2mzbhiTiPvtrqgTQJ+zqiAq2iPIPilpaQ==",
484+ "cpu": [
485+ "arm"
486+ ],
487+ "dev": true,
488+ "license": "MIT",
489+ "optional": true,
490+ "os": [
491+ "android"
492+ ]
493+ },
494+ "node_modules/@rollup/rollup-android-arm64": {
495+ "version": "4.63.1",
496+ "resolved": "https://registry.npmjs.org/@rollup/rollup-android-arm64/-/rollup-android-arm64-4.63.1.tgz",
497+ "integrity": "sha512-cQ4nFQABN5cDvDpbvJ7bMStCpnaVxynZrRMfUJYgxcIk9Sh54FIO1vtfkg0B69REjER77ioZ/ov+eAApx/KmLQ==",
498+ "cpu": [
499+ "arm64"
500+ ],
501+ "dev": true,
502+ "license": "MIT",
503+ "optional": true,
504+ "os": [
505+ "android"
506+ ]
507+ },
508+ "node_modules/@rollup/rollup-darwin-arm64": {
509+ "version": "4.63.1",
510+ "resolved": "https://registry.npmjs.org/@rollup/rollup-darwin-arm64/-/rollup-darwin-arm64-4.63.1.tgz",
511+ "integrity": "sha512-FQNqd1lRy/0QhDk3xeRIkSBiCpXCiDnZO3YLVdcDKN1UBiKToNftCzcXYNLshmPDUMlu2TdeS8tGcsU6f3YF1Q==",
512+ "cpu": [
513+ "arm64"
514+ ],
515+ "dev": true,
516+ "license": "MIT",
517+ "optional": true,
518+ "os": [
519+ "darwin"
520+ ]
521+ },
522+ "node_modules/@rollup/rollup-darwin-x64": {
523+ "version": "4.63.1",
524+ "resolved": "https://registry.npmjs.org/@rollup/rollup-darwin-x64/-/rollup-darwin-x64-4.63.1.tgz",
525+ "integrity": "sha512-pvD16V939D3CloK0+qikpGaxiPrDUXTe7Y5cWOMkMSy7m1cawa8EGy/kXYi/G/cKAC4HDAbSnzCIk1WmsoOKXg==",
526+ "cpu": [
527+ "x64"
528+ ],
529+ "dev": true,
530+ "license": "MIT",
531+ "optional": true,
532+ "os": [
533+ "darwin"
534+ ]
535+ },
536+ "node_modules/@rollup/rollup-freebsd-arm64": {
537+ "version": "4.63.1",
538+ "resolved": "https://registry.npmjs.org/@rollup/rollup-freebsd-arm64/-/rollup-freebsd-arm64-4.63.1.tgz",
539+ "integrity": "sha512-pcFGeL2345VwdTnJhA6zLbew+YgWB0qBG2+dMtXjCicf6+rm6kO6cOoh5VnTe0ZMrMRgRyuHmCJxZWrIdzYuOw==",
540+ "cpu": [
541+ "arm64"
542+ ],
543+ "dev": true,
544+ "license": "MIT",
545+ "optional": true,
546+ "os": [
547+ "freebsd"
548+ ]
549+ },
550+ "node_modules/@rollup/rollup-freebsd-x64": {
551+ "version": "4.63.1",
552+ "resolved": "https://registry.npmjs.org/@rollup/rollup-freebsd-x64/-/rollup-freebsd-x64-4.63.1.tgz",
553+ "integrity": "sha512-mRJlqSRulVzcKq/LKA6ICSIc3K/l4fzlVn/gePn2nXIHy8seRi5z/eeRE0d/XMBxcMldiXtQTSpRj0tkkC3g8Q==",
554+ "cpu": [
555+ "x64"
556+ ],
557+ "dev": true,
558+ "license": "MIT",
559+ "optional": true,
560+ "os": [
561+ "freebsd"
562+ ]
563+ },
564+ "node_modules/@rollup/rollup-linux-arm-gnueabihf": {
565+ "version": "4.63.1",
566+ "resolved": "https://registry.npmjs.org/@rollup/rollup-linux-arm-gnueabihf/-/rollup-linux-arm-gnueabihf-4.63.1.tgz",
567+ "integrity": "sha512-YDUNvVM85TI3g/1OpnqKP1h4NeW/j64DfWMf+G3M809xNk1bJSnpFp4sh83NpmVE5DXnkh8ULor4LTVZKoYLHw==",
568+ "cpu": [
569+ "arm"
570+ ],
571+ "dev": true,
572+ "libc": [
573+ "glibc"
574+ ],
575+ "license": "MIT",
576+ "optional": true,
577+ "os": [
578+ "linux"
579+ ]
580+ },
581+ "node_modules/@rollup/rollup-linux-arm-musleabihf": {
582+ "version": "4.63.1",
583+ "resolved": "https://registry.npmjs.org/@rollup/rollup-linux-arm-musleabihf/-/rollup-linux-arm-musleabihf-4.63.1.tgz",
584+ "integrity": "sha512-7Mcn71p9ZuQFAj+h+dhQXy/yeLePRS2yKRnmW1DijA9thKO5qap0GNOIQK4yQ6iP3SU0Mrb/yWo8h8vgRba8lw==",
585+ "cpu": [
586+ "arm"
587+ ],
588+ "dev": true,
589+ "libc": [
590+ "musl"
591+ ],
592+ "license": "MIT",
593+ "optional": true,
594+ "os": [
595+ "linux"
596+ ]
597+ },
598+ "node_modules/@rollup/rollup-linux-arm64-gnu": {
599+ "version": "4.63.1",
600+ "resolved": "https://registry.npmjs.org/@rollup/rollup-linux-arm64-gnu/-/rollup-linux-arm64-gnu-4.63.1.tgz",
601+ "integrity": "sha512-4YiLQTX6U4CSl0L9cluep9A9W6UmTfqBDc2/CH6wlu54pl4E7Jn3cOD8oxzvBDEGk/JMKgJ47C8g+radF7mwvg==",
602+ "cpu": [
603+ "arm64"
604+ ],
605+ "dev": true,
606+ "libc": [
607+ "glibc"
608+ ],
609+ "license": "MIT",
610+ "optional": true,
611+ "os": [
612+ "linux"
613+ ]
614+ },
615+ "node_modules/@rollup/rollup-linux-arm64-musl": {
616+ "version": "4.63.1",
617+ "resolved": "https://registry.npmjs.org/@rollup/rollup-linux-arm64-musl/-/rollup-linux-arm64-musl-4.63.1.tgz",
618+ "integrity": "sha512-2ra8F7w8OquwZN9z2/fKFnli69wa8PLwaVzRMIPGb13ByMJwC28Fbp8YcVGoUhlYMTt7j5j9bNgpysrN2UM+vw==",
619+ "cpu": [
620+ "arm64"
621+ ],
622+ "dev": true,
623+ "libc": [
624+ "musl"
625+ ],
626+ "license": "MIT",
627+ "optional": true,
628+ "os": [
629+ "linux"
630+ ]
631+ },
632+ "node_modules/@rollup/rollup-linux-loong64-gnu": {
633+ "version": "4.63.1",
634+ "resolved": "https://registry.npmjs.org/@rollup/rollup-linux-loong64-gnu/-/rollup-linux-loong64-gnu-4.63.1.tgz",
635+ "integrity": "sha512-Sy20ncyhjmBP0Ml+UvQbimjlk6VFgjW5uNP+qqwHB00mTE8Bl2C1TuHTlRwK2YoXeZbee5lP2XevBWVkAQAtSQ==",
636+ "cpu": [
637+ "loong64"
638+ ],
639+ "dev": true,
640+ "libc": [
641+ "glibc"
642+ ],
643+ "license": "MIT",
644+ "optional": true,
645+ "os": [
646+ "linux"
647+ ]
648+ },
649+ "node_modules/@rollup/rollup-linux-loong64-musl": {
650+ "version": "4.63.1",
651+ "resolved": "https://registry.npmjs.org/@rollup/rollup-linux-loong64-musl/-/rollup-linux-loong64-musl-4.63.1.tgz",
652+ "integrity": "sha512-noITLp8oNjYliPnGWmLyelIHwULGqbHloQHGw1rtxbWhTuWooRpnZarZQJ1y9EUC4szuCusCc+HEpUtxpIwYvA==",
653+ "cpu": [
654+ "loong64"
655+ ],
656+ "dev": true,
657+ "libc": [
658+ "musl"
659+ ],
660+ "license": "MIT",
661+ "optional": true,
662+ "os": [
663+ "linux"
664+ ]
665+ },
666+ "node_modules/@rollup/rollup-linux-ppc64-gnu": {
667+ "version": "4.63.1",
668+ "resolved": "https://registry.npmjs.org/@rollup/rollup-linux-ppc64-gnu/-/rollup-linux-ppc64-gnu-4.63.1.tgz",
669+ "integrity": "sha512-hlxxXd+F1mWiAcaFR7Sv9ZQT6m6UfI8+Vy/kFJzztq2pDMU/0wZ9sish0iszNZvsQDo8Gc0i5yuFEOz5dDf6fA==",
670+ "cpu": [
671+ "ppc64"
672+ ],
673+ "dev": true,
674+ "libc": [
675+ "glibc"
676+ ],
677+ "license": "MIT",
678+ "optional": true,
679+ "os": [
680+ "linux"
681+ ]
682+ },
683+ "node_modules/@rollup/rollup-linux-ppc64-musl": {
684+ "version": "4.63.1",
685+ "resolved": "https://registry.npmjs.org/@rollup/rollup-linux-ppc64-musl/-/rollup-linux-ppc64-musl-4.63.1.tgz",
686+ "integrity": "sha512-EF7OpqQTQ/BvGqLzUi4rEHuagCV9MugAUXSHemwPW5vxZ75RR+jxO/2j95Ph2dalMpFHSVECjRoioHZgA9zOYA==",
687+ "cpu": [
688+ "ppc64"
689+ ],
690+ "dev": true,
691+ "libc": [
692+ "musl"
693+ ],
694+ "license": "MIT",
695+ "optional": true,
696+ "os": [
697+ "linux"
698+ ]
699+ },
700+ "node_modules/@rollup/rollup-linux-riscv64-gnu": {
701+ "version": "4.63.1",
702+ "resolved": "https://registry.npmjs.org/@rollup/rollup-linux-riscv64-gnu/-/rollup-linux-riscv64-gnu-4.63.1.tgz",
703+ "integrity": "sha512-wQO3JesW9PRkwlabQ27y7sPfVOOTLRG73I4F2UYHG5PXun3J9U3y+b7ezVKSYbsvSKGQ1k1cq8Qlun4C9kLt3w==",
704+ "cpu": [
705+ "riscv64"
706+ ],
707+ "dev": true,
708+ "libc": [
709+ "glibc"
710+ ],
711+ "license": "MIT",
712+ "optional": true,
713+ "os": [
714+ "linux"
715+ ]
716+ },
717+ "node_modules/@rollup/rollup-linux-riscv64-musl": {
718+ "version": "4.63.1",
719+ "resolved": "https://registry.npmjs.org/@rollup/rollup-linux-riscv64-musl/-/rollup-linux-riscv64-musl-4.63.1.tgz",
720+ "integrity": "sha512-ouAGwhO6wHRXdnOVCOsB0tRFkA7nhNB2Nwax6oECXN0YiN8EYUTBAOudADOB1PI+yDL61TeNx/u7MVCzksNbkQ==",
721+ "cpu": [
722+ "riscv64"
723+ ],
724+ "dev": true,
725+ "libc": [
726+ "musl"
727+ ],
728+ "license": "MIT",
729+ "optional": true,
730+ "os": [
731+ "linux"
732+ ]
733+ },
734+ "node_modules/@rollup/rollup-linux-s390x-gnu": {
735+ "version": "4.63.1",
736+ "resolved": "https://registry.npmjs.org/@rollup/rollup-linux-s390x-gnu/-/rollup-linux-s390x-gnu-4.63.1.tgz",
737+ "integrity": "sha512-q2R38Sn+1J8RxhfJ+T54wSWmyKXWec+9jgDfqO2AtArEqHO5R2aeayp5H5OYLr5UYDVGsVaZPEFUooMhYCdz5A==",
738+ "cpu": [
739+ "s390x"
740+ ],
741+ "dev": true,
742+ "libc": [
743+ "glibc"
744+ ],
745+ "license": "MIT",
746+ "optional": true,
747+ "os": [
748+ "linux"
749+ ]
750+ },
751+ "node_modules/@rollup/rollup-linux-x64-gnu": {
752+ "version": "4.63.1",
753+ "resolved": "https://registry.npmjs.org/@rollup/rollup-linux-x64-gnu/-/rollup-linux-x64-gnu-4.63.1.tgz",
754+ "integrity": "sha512-gfI5T24WLLuFfSKw7Go/zDXjAAV0fny0swTaDv+WjK7vqcw4cRhFfdsyKL1n+ukI+ooBxn3bVQnyrn06WpI50w==",
755+ "cpu": [
756+ "x64"
757+ ],
758+ "dev": true,
759+ "libc": [
760+ "glibc"
761+ ],
762+ "license": "MIT",
763+ "optional": true,
764+ "os": [
765+ "linux"
766+ ]
767+ },
768+ "node_modules/@rollup/rollup-linux-x64-musl": {
769+ "version": "4.63.1",
770+ "resolved": "https://registry.npmjs.org/@rollup/rollup-linux-x64-musl/-/rollup-linux-x64-musl-4.63.1.tgz",
771+ "integrity": "sha512-4h6XqthmB4Hspji84wvgk+ElodTsGj+dbZqHJHHtKxj4mYq0ANSEEPX9ys3moJueqsRjwpaJYH7874Itwnj2ow==",
772+ "cpu": [
773+ "x64"
774+ ],
775+ "dev": true,
776+ "libc": [
777+ "musl"
778+ ],
779+ "license": "MIT",
780+ "optional": true,
781+ "os": [
782+ "linux"
783+ ]
784+ },
785+ "node_modules/@rollup/rollup-openbsd-x64": {
786+ "version": "4.63.1",
787+ "resolved": "https://registry.npmjs.org/@rollup/rollup-openbsd-x64/-/rollup-openbsd-x64-4.63.1.tgz",
788+ "integrity": "sha512-dlfCOa87o1VAYegLQ9EKilx2JCeRofiyPGhTCmqnuXZ6bMPiycO1rq1+sKoulAp7pGLIsTIw+1x5R+zgh5LhhA==",
789+ "cpu": [
790+ "x64"
791+ ],
792+ "dev": true,
793+ "license": "MIT",
794+ "optional": true,
795+ "os": [
796+ "openbsd"
797+ ]
798+ },
799+ "node_modules/@rollup/rollup-openharmony-arm64": {
800+ "version": "4.63.1",
801+ "resolved": "https://registry.npmjs.org/@rollup/rollup-openharmony-arm64/-/rollup-openharmony-arm64-4.63.1.tgz",
802+ "integrity": "sha512-cjkLbOlfcm3QGhMM1J5zaZjsw1GggbN6rw9UTSSRrPrR1KkcXnN7Uq9rPw34xImQ9VOY9GN+6u2Zj80B9ptkcw==",
803+ "cpu": [
804+ "arm64"
805+ ],
806+ "dev": true,
807+ "license": "MIT",
808+ "optional": true,
809+ "os": [
810+ "openharmony"
811+ ]
812+ },
813+ "node_modules/@rollup/rollup-win32-arm64-msvc": {
814+ "version": "4.63.1",
815+ "resolved": "https://registry.npmjs.org/@rollup/rollup-win32-arm64-msvc/-/rollup-win32-arm64-msvc-4.63.1.tgz",
816+ "integrity": "sha512-Li1KdUnWGE4N3e1F/B4RTB1ms+nG4WBgjByO46pkeBVX/2UBsY53xf5vK9WygVmnH3RwncIST7lkSdLSY6P9lg==",
817+ "cpu": [
818+ "arm64"
819+ ],
820+ "dev": true,
821+ "license": "MIT",
822+ "optional": true,
823+ "os": [
824+ "win32"
825+ ]
826+ },
827+ "node_modules/@rollup/rollup-win32-ia32-msvc": {
828+ "version": "4.63.1",
829+ "resolved": "https://registry.npmjs.org/@rollup/rollup-win32-ia32-msvc/-/rollup-win32-ia32-msvc-4.63.1.tgz",
830+ "integrity": "sha512-t4ZYOSoLTgwhuFMrmTMLx/+i1DQVK7HYqMc6kY46EApwi8X0nIVphzdNoThU3xt6n+N5urG1/gxBdCaKDLavfg==",
831+ "cpu": [
832+ "ia32"
833+ ],
834+ "dev": true,
835+ "license": "MIT",
836+ "optional": true,
837+ "os": [
838+ "win32"
839+ ]
840+ },
841+ "node_modules/@rollup/rollup-win32-x64-gnu": {
842+ "version": "4.63.1",
843+ "resolved": "https://registry.npmjs.org/@rollup/rollup-win32-x64-gnu/-/rollup-win32-x64-gnu-4.63.1.tgz",
844+ "integrity": "sha512-RgroPfMmKlD1RzSDxvwgcPiy2HNQKoYV7OmwIXDsk73uKW5t6B/V8KIy27SMv/FNXFo/oSBtWc9J0X7t91ezZg==",
845+ "cpu": [
846+ "x64"
847+ ],
848+ "dev": true,
849+ "license": "MIT",
850+ "optional": true,
851+ "os": [
852+ "win32"
853+ ]
854+ },
855+ "node_modules/@rollup/rollup-win32-x64-msvc": {
856+ "version": "4.63.1",
857+ "resolved": "https://registry.npmjs.org/@rollup/rollup-win32-x64-msvc/-/rollup-win32-x64-msvc-4.63.1.tgz",
858+ "integrity": "sha512-at8QVep6S3h5Y6gSbdGU06bRY5WJkf6WUduM9YtvYMbYhB1MOFfUgc6kehitQXzOtMSaT70q7f9ydPhpqu821w==",
859+ "cpu": [
860+ "x64"
861+ ],
862+ "dev": true,
863+ "license": "MIT",
864+ "optional": true,
865+ "os": [
866+ "win32"
867+ ]
868+ },
869+ "node_modules/@types/estree": {
870+ "version": "1.0.9",
871+ "resolved": "https://registry.npmjs.org/@types/estree/-/estree-1.0.9.tgz",
872+ "integrity": "sha512-GhdPgy1el4/ImP05X05Uw4cw2/M93BCUmnEvWZNStlCzEKME4Fkk+YpoA5OiHNQmoS7Cafb8Xa3Pya8m1Qrzeg==",
873+ "dev": true,
874+ "license": "MIT"
875+ },
876+ "node_modules/@types/node": {
877+ "version": "24.13.4",
878+ "resolved": "https://registry.npmjs.org/@types/node/-/node-24.13.4.tgz",
879+ "integrity": "sha512-YJ7EqCstVTzIr0fMr7qul/977en+pQHrfmuKIo6Zr9i75Be21dr3MovcfvGtyvi2HAUrRerWps5sMO9I7WaxDw==",
880+ "dev": true,
881+ "license": "MIT",
882+ "dependencies": {
883+ "undici-types": "~7.18.0"
884+ }
885+ },
886+ "node_modules/@webgpu/types": {
887+ "version": "0.1.72",
888+ "resolved": "https://registry.npmjs.org/@webgpu/types/-/types-0.1.72.tgz",
889+ "integrity": "sha512-0cF7RFM2edNoiIS1ODJp0/Gzv4/xSXhwoR0YCza+OWpJWtn4wmo9DvK91aLlH9+uUnwIriP7ZiC3WitmyhuzBw==",
890+ "dev": true,
891+ "license": "BSD-3-Clause"
892+ },
893+ "node_modules/esbuild": {
894+ "version": "0.25.12",
895+ "resolved": "https://registry.npmjs.org/esbuild/-/esbuild-0.25.12.tgz",
896+ "integrity": "sha512-bbPBYYrtZbkt6Os6FiTLCTFxvq4tt3JKall1vRwshA3fdVztsLAatFaZobhkBC8/BrPetoa0oksYoKXoG4ryJg==",
897+ "dev": true,
898+ "hasInstallScript": true,
899+ "license": "MIT",
900+ "bin": {
901+ "esbuild": "bin/esbuild"
902+ },
903+ "engines": {
904+ "node": ">=18"
905+ },
906+ "optionalDependencies": {
907+ "@esbuild/aix-ppc64": "0.25.12",
908+ "@esbuild/android-arm": "0.25.12",
909+ "@esbuild/android-arm64": "0.25.12",
910+ "@esbuild/android-x64": "0.25.12",
911+ "@esbuild/darwin-arm64": "0.25.12",
912+ "@esbuild/darwin-x64": "0.25.12",
913+ "@esbuild/freebsd-arm64": "0.25.12",
914+ "@esbuild/freebsd-x64": "0.25.12",
915+ "@esbuild/linux-arm": "0.25.12",
916+ "@esbuild/linux-arm64": "0.25.12",
917+ "@esbuild/linux-ia32": "0.25.12",
918+ "@esbuild/linux-loong64": "0.25.12",
919+ "@esbuild/linux-mips64el": "0.25.12",
920+ "@esbuild/linux-ppc64": "0.25.12",
921+ "@esbuild/linux-riscv64": "0.25.12",
922+ "@esbuild/linux-s390x": "0.25.12",
923+ "@esbuild/linux-x64": "0.25.12",
924+ "@esbuild/netbsd-arm64": "0.25.12",
925+ "@esbuild/netbsd-x64": "0.25.12",
926+ "@esbuild/openbsd-arm64": "0.25.12",
927+ "@esbuild/openbsd-x64": "0.25.12",
928+ "@esbuild/openharmony-arm64": "0.25.12",
929+ "@esbuild/sunos-x64": "0.25.12",
930+ "@esbuild/win32-arm64": "0.25.12",
931+ "@esbuild/win32-ia32": "0.25.12",
932+ "@esbuild/win32-x64": "0.25.12"
933+ }
934+ },
935+ "node_modules/fdir": {
936+ "version": "6.5.0",
937+ "resolved": "https://registry.npmjs.org/fdir/-/fdir-6.5.0.tgz",
938+ "integrity": "sha512-tIbYtZbucOs0BRGqPJkshJUYdL+SDH7dVM8gjy+ERp3WAUjLEFJE+02kanyHtwjWOnwrKYBiwAmM0p4kLJAnXg==",
939+ "dev": true,
940+ "license": "MIT",
941+ "engines": {
942+ "node": ">=12.0.0"
943+ },
944+ "peerDependencies": {
945+ "picomatch": "^3 || ^4"
946+ },
947+ "peerDependenciesMeta": {
948+ "picomatch": {
949+ "optional": true
950+ }
951+ }
952+ },
953+ "node_modules/fsevents": {
954+ "version": "2.3.3",
955+ "resolved": "https://registry.npmjs.org/fsevents/-/fsevents-2.3.3.tgz",
956+ "integrity": "sha512-5xoDfX+fL7faATnagmWPpbFtwh/R77WmMMqqHGS65C3vvB0YHrgF+B1YmZ3441tMj5n63k0212XNoJwzlhffQw==",
957+ "dev": true,
958+ "hasInstallScript": true,
959+ "license": "MIT",
960+ "optional": true,
961+ "os": [
962+ "darwin"
963+ ],
964+ "engines": {
965+ "node": "^8.16.0 || ^10.6.0 || >=11.0.0"
966+ }
967+ },
968+ "node_modules/nanoid": {
969+ "version": "3.3.18",
970+ "resolved": "https://registry.npmjs.org/nanoid/-/nanoid-3.3.18.tgz",
971+ "integrity": "sha512-DTg4MJbGMWkfi6VZFdNt2/caMbQy4Ou+Op/hJQvGEWcnVfoA1QA+xzRKAzw9jD6+GVOOeYr/mIcuDSdug6F6+w==",
972+ "dev": true,
973+ "funding": [
974+ {
975+ "type": "github",
976+ "url": "https://github.com/sponsors/ai"
977+ }
978+ ],
979+ "license": "MIT",
980+ "bin": {
981+ "nanoid": "bin/nanoid.cjs"
982+ },
983+ "engines": {
984+ "node": "^10 || ^12 || ^13.7 || ^14 || >=15.0.1"
985+ }
986+ },
987+ "node_modules/picocolors": {
988+ "version": "1.1.1",
989+ "resolved": "https://registry.npmjs.org/picocolors/-/picocolors-1.1.1.tgz",
990+ "integrity": "sha512-xceH2snhtb5M9liqDsmEw56le376mTZkEX/jEb/RxNFyegNul7eNslCXP9FDj/Lcu0X8KEyMceP2ntpaHrDEVA==",
991+ "dev": true,
992+ "license": "ISC"
993+ },
994+ "node_modules/picomatch": {
995+ "version": "4.0.7",
996+ "resolved": "https://registry.npmjs.org/picomatch/-/picomatch-4.0.7.tgz",
997+ "integrity": "sha512-qcJu88Q2IWqJsDD529JKMdwGm/dvInW4HvQnRwiH9JtihJvzGOscDtHE3x1pBKeUOTysQ8kVmLnJ2kJu7yhcGA==",
998+ "dev": true,
999+ "license": "MIT",
1000+ "engines": {
1001+ "node": ">=12"
1002+ },
1003+ "funding": {
1004+ "url": "https://github.com/sponsors/jonschlinkert"
1005+ }
1006+ },
1007+ "node_modules/playwright": {
1008+ "version": "1.63.0",
1009+ "resolved": "https://registry.npmjs.org/playwright/-/playwright-1.63.0.tgz",
1010+ "integrity": "sha512-+7ziBLidS4NaNCdt57SUDT+wYmmd5fmiQejUic/kb+YsYSCPyOOE9sebzMjNmQrsnNpDJqd4WHvV/8lfKfUDUg==",
1011+ "dev": true,
1012+ "license": "Apache-2.0",
1013+ "dependencies": {
1014+ "playwright-core": "1.63.0"
1015+ },
1016+ "bin": {
1017+ "playwright": "cli.js"
1018+ },
1019+ "engines": {
1020+ "node": ">=20"
1021+ }
1022+ },
1023+ "node_modules/playwright-core": {
1024+ "version": "1.63.0",
1025+ "resolved": "https://registry.npmjs.org/playwright-core/-/playwright-core-1.63.0.tgz",
1026+ "integrity": "sha512-rYCsBF/M5HjUch52bbtVONEFjv6Xu8sm8h72dNlR5bzIE1fvC/bxgspzkjSfU+MweEMmPM8KJebG6nnyxo5mCg==",
1027+ "dev": true,
1028+ "license": "Apache-2.0",
1029+ "bin": {
1030+ "playwright-core": "cli.js"
1031+ },
1032+ "engines": {
1033+ "node": ">=20"
1034+ }
1035+ },
1036+ "node_modules/postcss": {
1037+ "version": "8.5.28",
1038+ "resolved": "https://registry.npmjs.org/postcss/-/postcss-8.5.28.tgz",
1039+ "integrity": "sha512-RRuzqDtt5Y9h3quz5hWhK+TPnsmVs6WwSU6LkJMeY4HstUEDuYTG8UJSdawMRzmzAtV+KEoG8N3Qg2qLy5vM/A==",
1040+ "dev": true,
1041+ "funding": [
1042+ {
1043+ "type": "opencollective",
1044+ "url": "https://opencollective.com/postcss/"
1045+ },
1046+ {
1047+ "type": "tidelift",
1048+ "url": "https://tidelift.com/funding/github/npm/postcss"
1049+ },
1050+ {
1051+ "type": "github",
1052+ "url": "https://github.com/sponsors/ai"
1053+ }
1054+ ],
1055+ "license": "MIT",
1056+ "dependencies": {
1057+ "nanoid": "^3.3.18",
1058+ "picocolors": "^1.1.1",
1059+ "source-map-js": "^1.2.1"
1060+ },
1061+ "engines": {
1062+ "node": "^10 || ^12 || >=14"
1063+ }
1064+ },
1065+ "node_modules/rollup": {
1066+ "version": "4.63.1",
1067+ "resolved": "https://registry.npmjs.org/rollup/-/rollup-4.63.1.tgz",
1068+ "integrity": "sha512-3Df9jsstwhccuEfmAMi9l8XUh/GOkVObmFTU7CCVBysEbcOZLl84jCtaAZMcPiMz2EGKsATzQcU+Xr3n/wU6cg==",
1069+ "dev": true,
1070+ "license": "MIT",
1071+ "dependencies": {
1072+ "@types/estree": "1.0.9"
1073+ },
1074+ "bin": {
1075+ "rollup": "dist/bin/rollup"
1076+ },
1077+ "engines": {
1078+ "node": ">=18.0.0",
1079+ "npm": ">=8.0.0"
1080+ },
1081+ "optionalDependencies": {
1082+ "@napi-rs/lzma-linux-x64-gnu": "1.5.1",
1083+ "@rollup/rollup-android-arm-eabi": "4.63.1",
1084+ "@rollup/rollup-android-arm64": "4.63.1",
1085+ "@rollup/rollup-darwin-arm64": "4.63.1",
1086+ "@rollup/rollup-darwin-x64": "4.63.1",
1087+ "@rollup/rollup-freebsd-arm64": "4.63.1",
1088+ "@rollup/rollup-freebsd-x64": "4.63.1",
1089+ "@rollup/rollup-linux-arm-gnueabihf": "4.63.1",
1090+ "@rollup/rollup-linux-arm-musleabihf": "4.63.1",
1091+ "@rollup/rollup-linux-arm64-gnu": "4.63.1",
1092+ "@rollup/rollup-linux-arm64-musl": "4.63.1",
1093+ "@rollup/rollup-linux-loong64-gnu": "4.63.1",
1094+ "@rollup/rollup-linux-loong64-musl": "4.63.1",
1095+ "@rollup/rollup-linux-ppc64-gnu": "4.63.1",
1096+ "@rollup/rollup-linux-ppc64-musl": "4.63.1",
1097+ "@rollup/rollup-linux-riscv64-gnu": "4.63.1",
1098+ "@rollup/rollup-linux-riscv64-musl": "4.63.1",
1099+ "@rollup/rollup-linux-s390x-gnu": "4.63.1",
1100+ "@rollup/rollup-linux-x64-gnu": "4.63.1",
1101+ "@rollup/rollup-linux-x64-musl": "4.63.1",
1102+ "@rollup/rollup-openbsd-x64": "4.63.1",
1103+ "@rollup/rollup-openharmony-arm64": "4.63.1",
1104+ "@rollup/rollup-win32-arm64-msvc": "4.63.1",
1105+ "@rollup/rollup-win32-ia32-msvc": "4.63.1",
1106+ "@rollup/rollup-win32-x64-gnu": "4.63.1",
1107+ "@rollup/rollup-win32-x64-msvc": "4.63.1",
1108+ "fsevents": "~2.3.2"
1109+ }
1110+ },
1111+ "node_modules/source-map-js": {
1112+ "version": "1.2.1",
1113+ "resolved": "https://registry.npmjs.org/source-map-js/-/source-map-js-1.2.1.tgz",
1114+ "integrity": "sha512-UXWMKhLOwVKb728IUtQPXxfYU+usdybtUrK/8uGE8CQMvrhOpwvzDBwj0QhSL7MQc7vIsISBG8VQ8+IDQxpfQA==",
1115+ "dev": true,
1116+ "license": "BSD-3-Clause",
1117+ "engines": {
1118+ "node": ">=0.10.0"
1119+ }
1120+ },
1121+ "node_modules/tinyglobby": {
1122+ "version": "0.2.17",
1123+ "resolved": "https://registry.npmjs.org/tinyglobby/-/tinyglobby-0.2.17.tgz",
1124+ "integrity": "sha512-wXR/dYpcqKmfWpEdZjiKJOwCNFndD0DMnrW/cYjVGttEkBfVgcLFHoNrlj47mjOVic9yyNu65alsgF4NQyTa2g==",
1125+ "dev": true,
1126+ "license": "MIT",
1127+ "dependencies": {
1128+ "fdir": "^6.5.0",
1129+ "picomatch": "^4.0.4"
1130+ },
1131+ "engines": {
1132+ "node": ">=12.0.0"
1133+ },
1134+ "funding": {
1135+ "url": "https://github.com/sponsors/SuperchupuDev"
1136+ }
1137+ },
1138+ "node_modules/typescript": {
1139+ "version": "5.9.3",
1140+ "resolved": "https://registry.npmjs.org/typescript/-/typescript-5.9.3.tgz",
1141+ "integrity": "sha512-jl1vZzPDinLr9eUt3J/t7V6FgNEw9QjvBPdysz9KfQDD41fQrC2Y4vKQdiaUpFT4bXlb1RHhLpp8wtm6M5TgSw==",
1142+ "dev": true,
1143+ "license": "Apache-2.0",
1144+ "bin": {
1145+ "tsc": "bin/tsc",
1146+ "tsserver": "bin/tsserver"
1147+ },
1148+ "engines": {
1149+ "node": ">=14.17"
1150+ }
1151+ },
1152+ "node_modules/undici-types": {
1153+ "version": "7.18.2",
1154+ "resolved": "https://registry.npmjs.org/undici-types/-/undici-types-7.18.2.tgz",
1155+ "integrity": "sha512-AsuCzffGHJybSaRrmr5eHr81mwJU3kjw6M+uprWvCXiNeN9SOGwQ3Jn8jb8m3Z6izVgknn1R0FTCEAP2QrLY/w==",
1156+ "dev": true,
1157+ "license": "MIT"
1158+ },
1159+ "node_modules/vite": {
1160+ "version": "6.4.3",
1161+ "resolved": "https://registry.npmjs.org/vite/-/vite-6.4.3.tgz",
1162+ "integrity": "sha512-NTKlcQjlAK7MlQoyb6LgaqHc8sso/pVyUJYWMws3jg21uTJw/LddqIFPcPqP6PzpgbIcZyKI85sFE4HBrQDA8A==",
1163+ "dev": true,
1164+ "license": "MIT",
1165+ "dependencies": {
1166+ "esbuild": "^0.25.0",
1167+ "fdir": "^6.4.4",
1168+ "picomatch": "^4.0.2",
1169+ "postcss": "^8.5.3",
1170+ "rollup": "^4.34.9",
1171+ "tinyglobby": "^0.2.13"
1172+ },
1173+ "bin": {
1174+ "vite": "bin/vite.js"
1175+ },
1176+ "engines": {
1177+ "node": "^18.0.0 || ^20.0.0 || >=22.0.0"
1178+ },
1179+ "funding": {
1180+ "url": "https://github.com/vitejs/vite?sponsor=1"
1181+ },
1182+ "optionalDependencies": {
1183+ "fsevents": "~2.3.3"
1184+ },
1185+ "peerDependencies": {
1186+ "@types/node": "^18.0.0 || ^20.0.0 || >=22.0.0",
1187+ "jiti": ">=1.21.0",
1188+ "less": "*",
1189+ "lightningcss": "^1.21.0",
1190+ "sass": "*",
1191+ "sass-embedded": "*",
1192+ "stylus": "*",
1193+ "sugarss": "*",
1194+ "terser": "^5.16.0",
1195+ "tsx": "^4.8.1",
1196+ "yaml": "^2.4.2"
1197+ },
1198+ "peerDependenciesMeta": {
1199+ "@types/node": {
1200+ "optional": true
1201+ },
1202+ "jiti": {
1203+ "optional": true
1204+ },
1205+ "less": {
1206+ "optional": true
1207+ },
1208+ "lightningcss": {
1209+ "optional": true
1210+ },
1211+ "sass": {
1212+ "optional": true
1213+ },
1214+ "sass-embedded": {
1215+ "optional": true
1216+ },
1217+ "stylus": {
1218+ "optional": true
1219+ },
1220+ "sugarss": {
1221+ "optional": true
1222+ },
1223+ "terser": {
1224+ "optional": true
1225+ },
1226+ "tsx": {
1227+ "optional": true
1228+ },
1229+ "yaml": {
1230+ "optional": true
1231+ }
1232+ }
1233+ }
1234+ }
1235+}
package.jsonadded+22−0View file
@@ -0,0 +1,22 @@
1+{
2+ "name": "steelsim31",
3+ "version": "0.1.0",
4+ "private": true,
5+ "type": "module",
6+ "description": "Browser simulation of the sound of a hanging stainless steel sheet: geometrically nonlinear plate (Foppl-von Karman) stepped explicitly with WebGPU.",
7+ "scripts": {
8+ "dev": "vite",
9+ "build": "tsc --noEmit && vite build",
10+ "preview": "vite preview",
11+ "test": "node --test tests/*.test.ts",
12+ "test:cpu": "node --test tests/cpuRef.test.ts",
13+ "test:gpu": "node --test tests/gpu.test.ts"
14+ },
15+ "devDependencies": {
16+ "@types/node": "^24.0.0",
17+ "@webgpu/types": "^0.1.60",
18+ "playwright": "^1.50.0",
19+ "typescript": "^5.6.0",
20+ "vite": "^6.0.0"
21+ }
22+}
src/audio/fft.tsadded+46−0View file
@@ -0,0 +1,46 @@
1+// In-place radix-2 complex FFT.
2+
3+export class FFT {
4+ private cos: Float32Array;
5+ private sin: Float32Array;
6+ private rev: Uint32Array;
7+ constructor(readonly n: number) {
8+ if (n & (n - 1)) throw new Error("FFT size must be a power of two");
9+ this.cos = new Float32Array(n / 2);
10+ this.sin = new Float32Array(n / 2);
11+ for (let k = 0; k < n / 2; k++) {
12+ this.cos[k] = Math.cos((2 * Math.PI * k) / n);
13+ this.sin[k] = -Math.sin((2 * Math.PI * k) / n);
14+ }
15+ this.rev = new Uint32Array(n);
16+ const bits = Math.log2(n);
17+ for (let i = 0; i < n; i++) {
18+ let r = 0;
19+ for (let b = 0; b < bits; b++) r |= ((i >> b) & 1) << (bits - 1 - b);
20+ this.rev[i] = r;
21+ }
22+ }
23+ transform(re: Float32Array, im: Float32Array): void {
24+ const n = this.n;
25+ for (let i = 0; i < n; i++) {
26+ const r = this.rev[i];
27+ if (r > i) {
28+ let t = re[i]; re[i] = re[r]; re[r] = t;
29+ t = im[i]; im[i] = im[r]; im[r] = t;
30+ }
31+ }
32+ for (let size = 2; size <= n; size <<= 1) {
33+ const half = size >> 1, step = n / size;
34+ for (let start = 0; start < n; start += size) {
35+ for (let k = 0; k < half; k++) {
36+ const c = this.cos[k * step], s = this.sin[k * step];
37+ const a = start + k, b = a + half;
38+ const tr = re[b] * c - im[b] * s;
39+ const ti = re[b] * s + im[b] * c;
40+ re[b] = re[a] - tr; im[b] = im[a] - ti;
41+ re[a] += tr; im[a] += ti;
42+ }
43+ }
44+ }
45+ }
46+}
src/audio/level.tsadded+43−0View file
@@ -0,0 +1,43 @@
1+// Playback level: absolute (fixed pascal-to-full-scale mapping with a user gain),
2+// loudness normalization, or peak normalization. A soft clipper bounds the output.
3+
4+export type LevelMode = "absolute" | "loudness" | "peak";
5+
6+/** Linear gain applied to the pressure signal [Pa] for the given mode. */
7+export function levelGain(samples: Float32Array, fs: number, mode: LevelMode, gainDb: number, peak = 0.95): number {
8+ if (samples.length === 0) return 0;
9+ if (mode === "absolute") {
10+ // 0 dBFS corresponds to 1 Pa at 0 dB gain.
11+ return 10 ** (gainDb / 20);
12+ }
13+ if (mode === "peak") {
14+ let max = 0;
15+ for (let k = 0; k < samples.length; k++) max = Math.max(max, Math.abs(samples[k]));
16+ return max > 0 ? peak / max : 0;
17+ }
18+ // Loudness: scale the RMS of the sustained part (after the first 50 ms) to -16 dBFS.
19+ const start = Math.min(samples.length - 1, Math.round(0.05 * fs));
20+ let sum = 0;
21+ for (let k = start; k < samples.length; k++) sum += samples[k] * samples[k];
22+ const rms = Math.sqrt(sum / Math.max(1, samples.length - start));
23+ return rms > 0 ? 10 ** (-16 / 20) / rms : 0;
24+}
25+
26+/** Peak level [dBFS] the signal would reach before clipping. */
27+export function peakDbfs(samples: Float32Array, gain: number): number {
28+ let max = 0;
29+ for (let k = 0; k < samples.length; k++) max = Math.max(max, Math.abs(samples[k]));
30+ return max > 0 ? 20 * Math.log10(max * gain) : -Infinity;
31+}
32+
33+/** Returns a copy of `samples` scaled for playback or export, soft-clipped to [-peak, peak]. */
34+export function normalize(samples: Float32Array, fs: number, mode: LevelMode, gainDb: number, peak = 0.95): Float32Array {
35+ const out = new Float32Array(samples.length);
36+ const g = levelGain(samples, fs, mode, gainDb, peak);
37+ if (mode === "peak") {
38+ for (let k = 0; k < samples.length; k++) out[k] = samples[k] * g;
39+ return out;
40+ }
41+ for (let k = 0; k < samples.length; k++) out[k] = peak * Math.tanh((samples[k] * g) / peak);
42+ return out;
43+}
src/audio/player.tsadded+54−0View file
@@ -0,0 +1,54 @@
1+// Web Audio playback of a Float32Array signal with position reporting.
2+
3+import { normalize, type LevelMode } from "./level.ts";
4+
5+export class Player {
6+ private ctx: AudioContext | null = null;
7+ private source: AudioBufferSourceNode | null = null;
8+ private startedAt = 0;
9+ private offset = 0;
10+ private length = 0;
11+ onEnded: (() => void) | null = null;
12+
13+ get playing(): boolean {
14+ return this.source !== null;
15+ }
16+
17+ /** Current playback position [s], or null when stopped. */
18+ position(): number | null {
19+ if (!this.source || !this.ctx) return null;
20+ return Math.min(this.length, this.offset + this.ctx.currentTime - this.startedAt);
21+ }
22+
23+ play(samples: Float32Array, fs: number, offset = 0, mode: LevelMode = "peak", gainDb = 0): void {
24+ this.stop();
25+ if (samples.length === 0) return;
26+ if (!this.ctx) this.ctx = new AudioContext({ sampleRate: fs });
27+ const ctx = this.ctx;
28+ const buffer = ctx.createBuffer(1, samples.length, fs);
29+ buffer.getChannelData(0).set(normalize(samples, fs, mode, gainDb, 0.9));
30+ const src = ctx.createBufferSource();
31+ src.buffer = buffer;
32+ src.connect(ctx.destination);
33+ src.onended = () => {
34+ if (this.source === src) {
35+ this.source = null;
36+ this.onEnded?.();
37+ }
38+ };
39+ this.length = samples.length / fs;
40+ this.offset = Math.max(0, Math.min(offset, this.length));
41+ this.startedAt = ctx.currentTime;
42+ this.source = src;
43+ src.start(0, this.offset);
44+ void ctx.resume();
45+ }
46+
47+ stop(): void {
48+ if (this.source) {
49+ const s = this.source;
50+ this.source = null;
51+ try { s.stop(); } catch { /* already stopped */ }
52+ }
53+ }
54+}
src/audio/spectrogram.tsadded+233−0View file
@@ -0,0 +1,233 @@
1+// Incremental STFT spectrogram drawn on a 2D canvas.
2+
3+import { FFT } from "./fft.ts";
4+
5+export interface SpectrogramOptions {
6+ windowSize: number;
7+ hop: number;
8+ fmax: number; // Hz
9+ dbFloor: number; // dB below the peak
10+ logFrequency: boolean;
11+}
12+
13+export class Spectrogram {
14+ private fft: FFT;
15+ private window: Float32Array;
16+ private re: Float32Array;
17+ private im: Float32Array;
18+ /** Magnitude (dB) columns, each of length windowSize / 2. */
19+ private columns: Float32Array[] = [];
20+ private consumed = 0;
21+ private peakDb = -Infinity;
22+ private opts: SpectrogramOptions;
23+
24+ constructor(opts: SpectrogramOptions) {
25+ this.opts = opts;
26+ this.fft = new FFT(opts.windowSize);
27+ this.window = new Float32Array(opts.windowSize);
28+ for (let k = 0; k < opts.windowSize; k++) this.window[k] = 0.5 - 0.5 * Math.cos((2 * Math.PI * k) / opts.windowSize);
29+ this.re = new Float32Array(opts.windowSize);
30+ this.im = new Float32Array(opts.windowSize);
31+ }
32+
33+ get options(): SpectrogramOptions {
34+ return this.opts;
35+ }
36+
37+ /** Changes options; recomputes from scratch when the STFT parameters change. */
38+ configure(opts: Partial<SpectrogramOptions>, audio: Float32Array, nSamples: number): void {
39+ const next = { ...this.opts, ...opts };
40+ const recompute = next.windowSize !== this.opts.windowSize || next.hop !== this.opts.hop;
41+ this.opts = next;
42+ if (recompute) {
43+ this.fft = new FFT(next.windowSize);
44+ this.window = new Float32Array(next.windowSize);
45+ for (let k = 0; k < next.windowSize; k++) this.window[k] = 0.5 - 0.5 * Math.cos((2 * Math.PI * k) / next.windowSize);
46+ this.re = new Float32Array(next.windowSize);
47+ this.im = new Float32Array(next.windowSize);
48+ this.reset();
49+ this.update(audio, nSamples);
50+ }
51+ }
52+
53+ reset(): void {
54+ this.columns = [];
55+ this.consumed = 0;
56+ this.peakDb = -Infinity;
57+ }
58+
59+ /** Adds columns for all complete windows in audio[0..nSamples). */
60+ update(audio: Float32Array, nSamples: number): void {
61+ const { windowSize, hop } = this.opts;
62+ while (this.consumed + windowSize <= nSamples) {
63+ for (let k = 0; k < windowSize; k++) {
64+ this.re[k] = audio[this.consumed + k] * this.window[k];
65+ this.im[k] = 0;
66+ }
67+ this.fft.transform(this.re, this.im);
68+ const col = new Float32Array(windowSize / 2);
69+ for (let k = 0; k < windowSize / 2; k++) {
70+ const mag2 = this.re[k] * this.re[k] + this.im[k] * this.im[k];
71+ const db = 10 * Math.log10(mag2 + 1e-30);
72+ col[k] = db;
73+ if (db > this.peakDb) this.peakDb = db;
74+ }
75+ this.columns.push(col);
76+ this.consumed += hop;
77+ }
78+ }
79+
80+ get columnCount(): number {
81+ return this.columns.length;
82+ }
83+
84+ private imageKey = "";
85+ private image: ImageData | null = null;
86+ private drawKey = "";
87+
88+ /** Forces the next draw to recompute everything. */
89+ invalidate(): void {
90+ this.imageKey = "";
91+ this.drawKey = "";
92+ }
93+
94+ /**
95+ * Draws the spectrogram; `duration` [s] fixes the time axis, `playhead` [s] draws a cursor.
96+ * The plot image is cached and only recomputed when the data, options, or size change,
97+ * and nothing is drawn at all when the previous frame is still current.
98+ */
99+ draw(canvas: HTMLCanvasElement, fs: number, duration: number, playhead: number | null): void {
100+ const dpr = window.devicePixelRatio || 1;
101+ const W = Math.max(1, Math.floor(canvas.clientWidth * dpr));
102+ const H = Math.max(1, Math.floor(canvas.clientHeight * dpr));
103+ const o = this.opts;
104+ const imageKey = `${W}|${H}|${this.columns.length}|${duration}|${fs}|${o.windowSize}|${o.hop}|${o.fmax}|${o.dbFloor}|${o.logFrequency}`;
105+ const drawKey = `${imageKey}|${playhead === null ? "-" : playhead.toFixed(4)}`;
106+ if (drawKey === this.drawKey) return;
107+ this.drawKey = drawKey;
108+ if (canvas.width !== W || canvas.height !== H) {
109+ canvas.width = W;
110+ canvas.height = H;
111+ }
112+ const ctx = canvas.getContext("2d")!;
113+ // Whole pixels throughout: a fractional device pixel ratio would otherwise
114+ // make the ImageData and the pixel buffer disagree in size.
115+ const left = Math.round(52 * dpr), bottom = Math.round(22 * dpr), top = Math.round(6 * dpr), right = Math.round(8 * dpr);
116+ const pw = Math.max(1, Math.floor(W - left - right)), ph = Math.max(1, Math.floor(H - top - bottom));
117+ ctx.fillStyle = "#1c1e22";
118+ ctx.fillRect(0, 0, W, H);
119+ const { fmax, logFrequency } = this.opts;
120+ const fmin = logFrequency ? 20 : 0;
121+ const logMin = Math.log(fmin || 1), logMax = Math.log(fmax);
122+ if (imageKey !== this.imageKey || !this.image) {
123+ this.image = this.renderImage(ctx, pw, ph, fs, duration);
124+ this.imageKey = imageKey;
125+ }
126+ ctx.putImageData(this.image, left, top);
127+ this.drawAxes(ctx, dpr, left, top, pw, ph, fmin, fmax, logMin, logMax, logFrequency, duration);
128+ if (playhead !== null && playhead >= 0 && playhead <= duration) {
129+ const x = left + (playhead / duration) * pw;
130+ ctx.fillStyle = "rgba(255,255,255,0.8)";
131+ ctx.fillRect(x, top, Math.max(1, dpr), ph);
132+ }
133+ }
134+
135+ private renderImage(ctx: CanvasRenderingContext2D, pw: number, ph: number, fs: number, duration: number): ImageData {
136+ const img = ctx.createImageData(pw, ph);
137+ img.data.set(spectrogramPixels(this.columns, this.peakDb, this.opts, fs, duration, pw, ph));
138+ return img;
139+ }
140+
141+ /** STFT columns (dB) and their peak, for renderers outside the browser. */
142+ get data(): { columns: Float32Array[]; peakDb: number } {
143+ return { columns: this.columns, peakDb: this.peakDb };
144+ }
145+
146+ private drawAxes(ctx: CanvasRenderingContext2D, dpr: number, left: number, top: number, pw: number, ph: number, fmin: number, fmax: number, logMin: number, logMax: number, logFrequency: boolean, duration: number): void {
147+ ctx.fillStyle = "#c8cbd0";
148+ ctx.font = `${11 * dpr}px system-ui, sans-serif`;
149+ ctx.textAlign = "right";
150+ ctx.textBaseline = "middle";
151+ const fticks = logFrequency ? [20, 50, 100, 200, 500, 1000, 2000, 5000, 10000, 20000].filter((f) => f >= fmin && f <= fmax) : niceTicks(0, fmax, 6);
152+ for (const f of fticks) {
153+ const frac = logFrequency ? (Math.log(f) - logMin) / (logMax - logMin) : (f - fmin) / (fmax - fmin);
154+ const y = top + (1 - frac) * ph;
155+ ctx.fillText(f >= 1000 ? `${f / 1000}k` : `${f}`, left - 4 * dpr, y);
156+ ctx.fillRect(left - 3 * dpr, y, 3 * dpr, 1);
157+ }
158+ ctx.textAlign = "center";
159+ ctx.textBaseline = "top";
160+ for (const t of niceTicks(0, duration, 8)) {
161+ const x = left + (t / duration) * pw;
162+ ctx.fillText(`${+t.toFixed(2)}`, x, top + ph + 4 * dpr);
163+ ctx.fillRect(x, top + ph, 1, 3 * dpr);
164+ }
165+ ctx.save();
166+ ctx.translate(12 * dpr, top + ph / 2);
167+ ctx.rotate(-Math.PI / 2);
168+ ctx.textBaseline = "middle";
169+ ctx.fillText("Hz", 0, 0);
170+ ctx.restore();
171+ }
172+}
173+
174+/** RGBA pixels (row-major, top row first) of the spectrogram plot area. */
175+export function spectrogramPixels(columns: Float32Array[], peakDb: number, opts: SpectrogramOptions, fs: number, duration: number, pw: number, ph: number, peakRef?: number): Uint8ClampedArray {
176+ const { windowSize, hop, fmax, dbFloor, logFrequency } = opts;
177+ const binHz = fs / windowSize;
178+ const fmin = logFrequency ? 20 : 0;
179+ const data = new Uint8ClampedArray(pw * ph * 4);
180+ const cols = columns.length;
181+ const colPerPx = (duration * fs) / hop / pw;
182+ const logMin = Math.log(fmin || 1), logMax = Math.log(fmax);
183+ const ref = peakRef ?? peakDb;
184+ for (let py = 0; py < ph; py++) {
185+ const frac = 1 - py / ph;
186+ const f = logFrequency ? Math.exp(logMin + frac * (logMax - logMin)) : fmin + frac * (fmax - fmin);
187+ const bin = Math.min(windowSize / 2 - 1, Math.max(0, Math.round(f / binHz)));
188+ for (let px = 0; px < pw; px++) {
189+ const c0 = Math.floor(px * colPerPx);
190+ const c1 = Math.max(c0 + 1, Math.floor((px + 1) * colPerPx));
191+ let v = -Infinity;
192+ for (let c = c0; c < c1 && c < cols; c++) v = Math.max(v, columns[c][bin]);
193+ const o = (py * pw + px) * 4;
194+ if (c0 >= cols) {
195+ data[o] = 28; data[o + 1] = 30; data[o + 2] = 34; data[o + 3] = 255;
196+ continue;
197+ }
198+ const t = ref > -200 ? Math.max(0, Math.min(1, (v - ref + dbFloor) / dbFloor)) : 0;
199+ const [r, g, b] = inferno(t);
200+ data[o] = r; data[o + 1] = g; data[o + 2] = b; data[o + 3] = 255;
201+ }
202+ }
203+ return data;
204+}
205+
206+function niceTicks(lo: number, hi: number, count: number): number[] {
207+ const span = hi - lo;
208+ if (span <= 0) return [lo];
209+ const raw = span / count;
210+ const mag = 10 ** Math.floor(Math.log10(raw));
211+ const norm = raw / mag;
212+ const step = (norm < 1.5 ? 1 : norm < 3.5 ? 2 : norm < 7.5 ? 5 : 10) * mag;
213+ const ticks: number[] = [];
214+ for (let t = Math.ceil(lo / step) * step; t <= hi + 1e-9; t += step) ticks.push(t);
215+ return ticks;
216+}
217+
218+// Compact polynomial approximation of the inferno colormap.
219+function inferno(t: number): [number, number, number] {
220+ const c0 = [0.0002189403691192265, 0.001651004631001012, -0.01948089843709184];
221+ const c1 = [0.1065134194856116, 0.5639564367884091, 3.932712388889277];
222+ const c2 = [11.60249308247187, -3.972853965665698, -15.9423940062914];
223+ const c3 = [-41.70399613139459, 17.43639888205313, 44.35414519872813];
224+ const c4 = [77.162935699427, -33.40235894210092, -81.80730925738993];
225+ const c5 = [-71.31942824499214, 32.62606426397723, 73.20951985803202];
226+ const c6 = [25.13112622477341, -12.24266895238567, -23.07032500287172];
227+ const out: number[] = [];
228+ for (let k = 0; k < 3; k++) {
229+ const v = c0[k] + t * (c1[k] + t * (c2[k] + t * (c3[k] + t * (c4[k] + t * (c5[k] + t * c6[k])))));
230+ out.push(Math.max(0, Math.min(255, Math.round(v * 255))));
231+ }
232+ return [out[0], out[1], out[2]];
233+}
src/audio/wav.tsadded+33−0View file
@@ -0,0 +1,33 @@
1+// 16-bit PCM mono WAV encoder.
2+
3+import { normalize, type LevelMode } from "./level.ts";
4+
5+/** Encodes the scaled signal as WAV bytes. */
6+export function encodeWavBytes(samples: Float32Array, fs: number, mode: LevelMode = "peak", gainDb = 0): Uint8Array<ArrayBuffer> {
7+ const scaled = normalize(samples, fs, mode, gainDb, 0.99);
8+ const buf = new ArrayBuffer(44 + samples.length * 2);
9+ const v = new DataView(buf);
10+ const str = (o: number, s: string) => { for (let k = 0; k < s.length; k++) v.setUint8(o + k, s.charCodeAt(k)); };
11+ str(0, "RIFF");
12+ v.setUint32(4, 36 + samples.length * 2, true);
13+ str(8, "WAVE");
14+ str(12, "fmt ");
15+ v.setUint32(16, 16, true);
16+ v.setUint16(20, 1, true);
17+ v.setUint16(22, 1, true);
18+ v.setUint32(24, fs, true);
19+ v.setUint32(28, fs * 2, true);
20+ v.setUint16(32, 2, true);
21+ v.setUint16(34, 16, true);
22+ str(36, "data");
23+ v.setUint32(40, samples.length * 2, true);
24+ for (let k = 0; k < samples.length; k++) {
25+ const s = Math.max(-1, Math.min(1, scaled[k]));
26+ v.setInt16(44 + 2 * k, Math.round(s * 32767), true);
27+ }
28+ return new Uint8Array(buf);
29+}
30+
31+export function encodeWav(samples: Float32Array, fs: number, mode: LevelMode = "peak", gainDb = 0): Blob {
32+ return new Blob([encodeWavBytes(samples, fs, mode, gainDb)], { type: "audio/wav" });
33+}
src/main.tsadded+489−0View file
@@ -0,0 +1,489 @@
1+// Application: parameter form, run loop, live and movie views, spectrogram,
2+// playback, WAV download, and the test hooks on window.steelsim.
3+
4+import { defaultParams, type PhysParams } from "./params.ts";
5+import { Simulation, type SimStatus } from "./sim/simulation.ts";
6+import { Renderer, type ViewMode } from "./render/renderer.ts";
7+import { Spectrogram } from "./audio/spectrogram.ts";
8+import { Player } from "./audio/player.ts";
9+import { encodeWav } from "./audio/wav.ts";
10+import { levelGain, peakDbfs, type LevelMode } from "./audio/level.ts";
11+import { ParamForm } from "./ui/form.ts";
12+import { runSelfTest, selfTestPassed, type SelfTestCase } from "./selftest.ts";
13+
14+const $ = <T extends HTMLElement>(id: string) => document.getElementById(id) as T;
15+
16+// ---------- State ----------
17+
18+const params: PhysParams = { ...defaultParams };
19+let form: ParamForm;
20+let device: GPUDevice | null = null;
21+let sim: Simulation | null = null;
22+let renderer: Renderer | null = null;
23+let running = false;
24+let paused = false;
25+let live = true;
26+let scrubFrame = 0;
27+let movieBuf: GPUBuffer | null = null;
28+let movieBufFrame = -1;
29+const spec = new Spectrogram({ windowSize: 2048, hop: 512, fmax: 5000, dbFloor: 80, logFrequency: false });
30+const player = new Player();
31+let lastStatsTime = 0, lastStatsSteps = 0, stepsPerSec = 0;
32+let duration = params.duration;
33+
34+// ---------- Controls ----------
35+
36+const btnStart = mkButton("Start", "primary");
37+const btnPause = mkButton("Pause");
38+const btnStop = mkButton("Stop");
39+const btnContinue = mkButton("Continue");
40+const durInput = document.createElement("input");
41+
42+function mkButton(text: string, cls = ""): HTMLButtonElement {
43+ const b = document.createElement("button");
44+ b.textContent = text;
45+ b.className = cls;
46+ return b;
47+}
48+
49+function buildControls(): void {
50+ const c = $("controls");
51+ const durRow = document.createElement("div");
52+ durRow.className = "field";
53+ durRow.style.width = "100%";
54+ const l = document.createElement("label");
55+ l.textContent = "Recorded duration";
56+ durInput.type = "number";
57+ durInput.step = "0.5";
58+ durInput.min = "0.1";
59+ durInput.value = String(duration);
60+ durInput.addEventListener("change", () => {
61+ const v = parseFloat(durInput.value);
62+ if (isFinite(v) && v > 0) {
63+ duration = v;
64+ params.duration = v;
65+ sim?.setDuration(v);
66+ updateButtons();
67+ }
68+ });
69+ const u = document.createElement("span");
70+ u.className = "unit";
71+ u.textContent = "s";
72+ durRow.append(l, durInput, u);
73+ c.append(durRow, btnStart, btnPause, btnStop, btnContinue);
74+ btnStart.addEventListener("click", () => void start());
75+ btnPause.addEventListener("click", () => {
76+ paused = !paused;
77+ btnPause.textContent = paused ? "Resume" : "Pause";
78+ if (!paused) void runLoop();
79+ });
80+ btnStop.addEventListener("click", () => {
81+ if (!sim) return;
82+ sim.stop();
83+ paused = false;
84+ btnPause.textContent = "Pause";
85+ if (!running) {
86+ form.setEnabled(true);
87+ updateButtons();
88+ }
89+ });
90+ btnContinue.addEventListener("click", () => {
91+ if (!sim) return;
92+ sim.setDuration(duration);
93+ if (sim.phase !== "done") void runLoop();
94+ updateButtons();
95+ });
96+ updateButtons();
97+}
98+
99+function updateButtons(): void {
100+ const s = sim?.status();
101+ btnPause.disabled = !sim || s?.phase === "done" || (!running && !paused);
102+ btnStop.disabled = !sim || s?.phase === "done";
103+ btnContinue.disabled = !sim || running || !sim.canContinue();
104+}
105+
106+// ---------- View toolbars ----------
107+
108+const viewMode = { value: "2d" as ViewMode };
109+const scrub = document.createElement("input");
110+const timeLabel = document.createElement("span");
111+const btnLive = mkButton("Live");
112+const btnPlay = mkButton("Play");
113+const btnWav = mkButton("Download .wav");
114+const levelSel = document.createElement("select");
115+const gainSlider = document.createElement("input");
116+const gainLabel = document.createElement("span");
117+const peakLabel = document.createElement("span");
118+let exagInput: HTMLInputElement;
119+let rangeInput: HTMLInputElement;
120+let shapeCheck: HTMLInputElement;
121+
122+function gainDb(): number {
123+ return parseFloat(gainSlider.value) || 0;
124+}
125+function levelMode(): LevelMode {
126+ return levelSel.value as LevelMode;
127+}
128+
129+function updateLevelReadout(): void {
130+ gainLabel.textContent = `${gainDb() >= 0 ? "+" : ""}${gainDb()} dB`;
131+ gainSlider.disabled = levelMode() !== "absolute";
132+ if (!sim || sim.samples === 0) { peakLabel.textContent = ""; return; }
133+ const audio = sim.recordedAudio();
134+ const db = peakDbfs(audio, levelGain(audio, sim.fs, levelMode(), gainDb()));
135+ peakLabel.textContent = isFinite(db) ? `peak ${db.toFixed(1)} dBFS${db > 0 ? " (clipped)" : ""}` : "";
136+ peakLabel.className = db > 0 ? "warn" : "sub";
137+ peakLabel.style.margin = "0";
138+}
139+
140+function buildViewToolbar(): void {
141+ const tb = $("viewToolbar");
142+ const seg = document.createElement("span");
143+ seg.className = "seg";
144+ const b2 = mkButton("2D", "on"), b3 = mkButton("3D");
145+ b2.addEventListener("click", () => { viewMode.value = "2d"; b2.classList.add("on"); b3.classList.remove("on"); });
146+ b3.addEventListener("click", () => { viewMode.value = "3d"; b3.classList.add("on"); b2.classList.remove("on"); });
147+ seg.append(b2, b3);
148+ exagInput = numInput(4, 0.5);
149+ rangeInput = numInput(0, 0.5);
150+ shapeCheck = document.createElement("input");
151+ shapeCheck.type = "checkbox";
152+ const hint = document.createElement("span");
153+ hint.className = "sub";
154+ hint.style.margin = "0";
155+ hint.textContent = "3D: drag to orbit, wheel to zoom";
156+ tb.append(seg, labeled("Exaggeration ×", exagInput), labeled("Color range (mm, 0 = auto)", rangeInput), labeled("3D: add initial shape", shapeCheck), hint);
157+}
158+
159+function buildMovieToolbar(): void {
160+ const tb = $("movieToolbar");
161+ scrub.type = "range";
162+ scrub.min = "0"; scrub.max = "0"; scrub.value = "0";
163+ timeLabel.className = "time";
164+ // Playback drives the movie; afterwards the view returns to where it was (live or scrubbed).
165+ let liveBeforePlay = true;
166+ const endPlayback = () => {
167+ btnPlay.textContent = "Play";
168+ if (liveBeforePlay) {
169+ live = true;
170+ updateLiveButton();
171+ }
172+ };
173+ scrub.addEventListener("input", () => {
174+ live = false;
175+ liveBeforePlay = false;
176+ scrubFrame = parseInt(scrub.value);
177+ if (player.playing) player.stop();
178+ updateLiveButton();
179+ });
180+ btnLive.addEventListener("click", () => { live = true; liveBeforePlay = true; updateLiveButton(); });
181+ btnPlay.addEventListener("click", () => {
182+ if (player.playing) { player.stop(); endPlayback(); return; }
183+ if (!sim || sim.samples === 0) return;
184+ liveBeforePlay = live;
185+ const offset = live ? 0 : scrubFrame / sim.fps;
186+ live = false;
187+ player.play(sim.recordedAudio(), sim.fs, offset, levelMode(), gainDb());
188+ btnPlay.textContent = "Stop";
189+ updateLiveButton();
190+ });
191+ player.onEnded = endPlayback;
192+ btnWav.addEventListener("click", () => {
193+ if (!sim || sim.samples === 0) return;
194+ const blob = encodeWav(sim.recordedAudio(), sim.fs, levelMode(), gainDb());
195+ const a = document.createElement("a");
196+ a.href = URL.createObjectURL(blob);
197+ a.download = `sheet_${sim.samples}samples_${sim.fs}Hz.wav`;
198+ a.click();
199+ setTimeout(() => URL.revokeObjectURL(a.href), 5000);
200+ });
201+ for (const [v, t] of [["absolute", "Absolute (0 dBFS = 1 Pa at 0 dB)"], ["loudness", "Loudness (RMS, soft clip)"], ["peak", "Peak"]]) {
202+ const o = document.createElement("option");
203+ o.value = v; o.textContent = t;
204+ levelSel.appendChild(o);
205+ }
206+ gainSlider.type = "range";
207+ gainSlider.min = "-40"; gainSlider.max = "60"; gainSlider.step = "1"; gainSlider.value = "20";
208+ gainSlider.style.width = "120px";
209+ gainSlider.style.flex = "none";
210+ const restart = () => {
211+ updateLevelReadout();
212+ const pos = player.position();
213+ if (pos !== null && sim) player.play(sim.recordedAudio(), sim.fs, pos, levelMode(), gainDb());
214+ };
215+ gainSlider.addEventListener("input", restart);
216+ levelSel.addEventListener("change", restart);
217+ tb.append(btnPlay, btnLive, scrub, timeLabel, labeled("Level", levelSel), labeled("Gain", gainSlider), gainLabel, peakLabel, btnWav);
218+ updateLevelReadout();
219+ updateLiveButton();
220+}
221+
222+function updateLiveButton(): void {
223+ btnLive.classList.toggle("on", live);
224+ btnLive.disabled = live;
225+}
226+
227+function buildSpecToolbar(): void {
228+ const tb = $("specToolbar");
229+ const win = document.createElement("select");
230+ for (const n of [512, 1024, 2048, 4096, 8192]) {
231+ const o = document.createElement("option");
232+ o.value = String(n); o.textContent = String(n);
233+ win.appendChild(o);
234+ }
235+ win.value = "2048";
236+ const hop = numInput(512, 64);
237+ const fmax = numInput(5000, 500);
238+ const floor = numInput(80, 5);
239+ const logf = document.createElement("input");
240+ logf.type = "checkbox";
241+ const apply = () => {
242+ if (!sim) return;
243+ spec.configure(
244+ { windowSize: parseInt(win.value), hop: Math.max(16, parseInt(hop.value)), fmax: Math.min(sim.fs / 2, parseFloat(fmax.value)), dbFloor: parseFloat(floor.value), logFrequency: logf.checked },
245+ sim.recordedAudio(), sim.samples,
246+ );
247+ };
248+ for (const el of [win, hop, fmax, floor, logf]) el.addEventListener("change", apply);
249+ tb.append(document.createTextNode("Spectrogram"), labeled("Window", win), labeled("Hop", hop), labeled("Max Hz", fmax), labeled("dB range", floor), labeled("Log freq.", logf));
250+}
251+
252+function numInput(v: number, step: number): HTMLInputElement {
253+ const i = document.createElement("input");
254+ i.type = "number"; i.value = String(v); i.step = String(step);
255+ return i;
256+}
257+function labeled(text: string, el: HTMLElement): HTMLLabelElement {
258+ const l = document.createElement("label");
259+ l.append(text, el);
260+ return l;
261+}
262+
263+// ---------- Simulation control ----------
264+
265+async function start(): Promise<void> {
266+ if (!device) return;
267+ player.stop();
268+ running = false;
269+ await new Promise((r) => setTimeout(r, 0));
270+ sim?.destroy();
271+ spec.reset();
272+ live = true;
273+ paused = false;
274+ btnPause.textContent = "Pause";
275+ params.duration = duration;
276+ sim = new Simulation(device, params);
277+ renderer = new Renderer(device, $("view"), sim.grid, params, sim.w0, sim.nw);
278+ movieBuf?.destroy();
279+ movieBuf = device.createBuffer({ size: Math.max(16, sim.grid.n * 4), usage: GPUBufferUsage.STORAGE | GPUBufferUsage.COPY_DST });
280+ movieBufFrame = -1;
281+ scrub.max = "0"; scrub.value = "0";
282+ lastStatsTime = performance.now(); lastStatsSteps = 0; stepsPerSec = 0;
283+ form.setEnabled(false);
284+ updateLiveButton();
285+ void runLoop();
286+}
287+
288+async function runLoop(): Promise<void> {
289+ if (!sim || running) return;
290+ running = true;
291+ updateButtons();
292+ try {
293+ while (sim.phase !== "done" && !paused) {
294+ await sim.runBatch();
295+ spec.update(sim.recordedAudio(), sim.samples);
296+ scrub.max = String(Math.max(0, sim.frames.length - 1));
297+ if (live) scrub.value = scrub.max;
298+ updateLevelReadout();
299+ const now = performance.now();
300+ if (now - lastStatsTime > 500) {
301+ stepsPerSec = ((sim.totalSteps - lastStatsSteps) * 1000) / (now - lastStatsTime);
302+ lastStatsTime = now; lastStatsSteps = sim.totalSteps;
303+ }
304+ }
305+ } catch (e) {
306+ console.error(e);
307+ $("info").innerHTML = `<div class="err">Simulation error: ${String(e)}</div>`;
308+ }
309+ running = false;
310+ if (sim.phase === "done") form.setEnabled(true);
311+ updateButtons();
312+}
313+
314+// ---------- Render loop ----------
315+
316+let frameErrorShown = false;
317+let lastRenderKey = "";
318+let lastInfoHtml = "";
319+
320+function frame(): void {
321+ requestAnimationFrame(frame);
322+ try {
323+ frameBody();
324+ } catch (e) {
325+ if (!frameErrorShown) { frameErrorShown = true; console.error("frame exception:", e); }
326+ }
327+}
328+
329+function frameBody(): void {
330+ if (!device || !sim || !renderer) return;
331+ renderer.mode = viewMode.value;
332+ renderer.exaggeration = parseFloat(exagInput.value) || 1;
333+ renderer.colorRange = (parseFloat(rangeInput.value) || 0) * 1e-3;
334+ renderer.showShape = shapeCheck.checked;
335+ const pos = player.position();
336+ if (pos !== null) {
337+ scrubFrame = Math.min(sim.frames.length - 1, Math.floor(pos * sim.fps));
338+ scrub.value = String(scrubFrame);
339+ }
340+ const st = sim.status();
341+ // Render only when something visible changed, so an idle or paused page does no work.
342+ const canvas = $<HTMLCanvasElement>("view");
343+ const cam = renderer.camera;
344+ const renderKey = [
345+ renderer.mode, renderer.exaggeration, renderer.colorRange, renderer.showShape,
346+ cam.yaw.toFixed(4), cam.pitch.toFixed(4), cam.dist.toFixed(4),
347+ live, scrubFrame, sim.frames.length, st.totalSteps, st.phase,
348+ canvas.clientWidth, canvas.clientHeight, window.devicePixelRatio,
349+ ].join("|");
350+ if (renderKey !== lastRenderKey) {
351+ lastRenderKey = renderKey;
352+ renderView(device, st);
353+ }
354+ spec.draw($("spec"), sim.fs, Math.max(duration, sim.samples / sim.fs), pos ?? (live ? null : scrubFrame / sim.fps));
355+ updateInfo(st);
356+}
357+
358+function renderView(dev: GPUDevice, st: SimStatus): void {
359+ if (!sim || !renderer) return;
360+ const enc = dev.createCommandEncoder();
361+ if (live || sim.frames.length === 0) {
362+ sim.encodeDisplay(enc);
363+ renderer.render(enc, sim.displayW, st.wmax);
364+ timeLabel.textContent = `${st.t >= 0 ? "" : "prep "}${st.t.toFixed(3)} s`;
365+ } else {
366+ const f = Math.max(0, Math.min(sim.frames.length - 1, scrubFrame));
367+ if (movieBufFrame !== f) {
368+ dev.queue.writeBuffer(movieBuf!, 0, sim.frames[f]);
369+ movieBufFrame = f;
370+ }
371+ let wmax = 0;
372+ const fr = sim.frames[f];
373+ for (let k = 0; k < fr.length; k += 7) wmax = Math.max(wmax, Math.abs(fr[k]));
374+ renderer.render(enc, movieBuf!, wmax);
375+ timeLabel.textContent = `${(f / sim.fps).toFixed(3)} s`;
376+ }
377+ dev.pushErrorScope("validation");
378+ dev.queue.submit([enc.finish()]);
379+ void dev.popErrorScope().then((e) => { if (e && !frameErrorShown) { frameErrorShown = true; console.error("frame validation error:", e.message); } });
380+}
381+
382+function updateInfo(st: SimStatus): void {
383+ if (!sim) return;
384+ const g = sim.grid;
385+ const { inside } = form.gridInfo();
386+ const spl = st.peakPressure > 0 ? 20 * Math.log10(st.peakPressure / 2e-5) : -Infinity;
387+ const phase = { push: "pushing", hold: "settling (held)", recording: "ringing (recording)", done: "done" }[st.phase];
388+ const rows: [string, string][] = [
389+ ["Phase", phase],
390+ ["Time", `${st.t.toFixed(3)} s`],
391+ ["Grid", `${g.nx} × ${g.ny} nodes${g.mask ? ` (${inside} in the sheet)` : ""}, ${(g.hx * 1e3).toFixed(2)} × ${(g.hy * 1e3).toFixed(2)} mm`],
392+ ["Time step", `${(sim.dt * 1e6).toFixed(3)} µs (${sim.stepsPerSample} per sample)`],
393+ ["Speed", `${(stepsPerSec / 1e3).toFixed(0)} k steps/s (${(stepsPerSec * sim.dt).toFixed(3)}× real time)`],
394+ ["Kinetic energy", `${st.ke.toExponential(2)} J`],
395+ ["Max |w|", `${(st.wmax * 1e3).toFixed(3)} mm`],
396+ ["Audio", `${(st.samples / sim.fs).toFixed(2)} s, ${st.frames} frames`],
397+ ["Peak pressure", st.peakPressure > 0 ? `${st.peakPressure.toFixed(3)} Pa (${spl.toFixed(1)} dB SPL)` : "–"],
398+ ];
399+ if (st.stick) {
400+ rows.push(["Stick", `${(st.stick.x * 1e3).toFixed(2)} mm, ${st.stick.v.toFixed(2)} m/s`]);
401+ rows.push(["Contact force", `${st.stick.force.toFixed(1)} N (max ${st.stick.maxForce.toFixed(1)} N)`]);
402+ }
403+ let html = rows.map(([a, b]) => `<div class="row"><span>${a}</span><span>${b}</span></div>`).join("");
404+ if (st.delayClamped) html += `<div class="warn">The listener is so close that some propagation delays were clamped.</div>`;
405+ if (!isFinite(st.wmax) || st.wmax > 1) html += `<div class="err">The solution has diverged. Increase the grid spacing (which lengthens the time step) or reduce the excitation.</div>`;
406+ if (html !== lastInfoHtml) {
407+ lastInfoHtml = html;
408+ $("info").innerHTML = html;
409+ }
410+}
411+
412+// ---------- Self-test ----------
413+
414+let selfTestResult: { done: boolean; passed?: boolean; cases?: SelfTestCase[]; error?: string } = { done: false };
415+
416+async function selfTest(): Promise<void> {
417+ if (!device) { selfTestResult = { done: true, passed: false, error: "no GPU device" }; return; }
418+ $("info").innerHTML = `<p class="sub">Running the GPU-versus-CPU self-test...</p>`;
419+ try {
420+ const cases = await runSelfTest(device);
421+ const passed = cases.every(selfTestPassed);
422+ selfTestResult = { done: true, passed, cases };
423+ const rows = cases.map((c) => `<div class="row"><span>${c.name}</span><span class="${selfTestPassed(c) ? "" : "err"}">acc ${c.accError.u.toExponential(1)} / ${c.accError.v.toExponential(1)} / ${c.accError.w.toExponential(1)}, pos ${c.positionError.toExponential(1)}</span></div>`);
424+ $("info").innerHTML = `<div class="${passed ? "row" : "err"}">Self-test ${passed ? "passed" : "FAILED"}</div>${rows.join("")}<p class="sub">Maximum relative errors of the u, v, w accelerations and of the compensated w update, GPU step versus float64 reference.</p>`;
425+ } catch (e) {
426+ selfTestResult = { done: true, passed: false, error: String(e) };
427+ $("info").innerHTML = `<div class="err">Self-test error: ${String(e)}</div>`;
428+ }
429+}
430+
431+// ---------- Init ----------
432+
433+async function init(): Promise<void> {
434+ form = new ParamForm($("params"), params);
435+ form.onChange = (key, isLive) => {
436+ if (!isLive || !sim) return;
437+ if (key === "sigma0" || key === "sigma1" || key === "sigmaIn") sim.setDamping(params.sigma0, params.sigma1, params.sigmaIn);
438+ if (key.startsWith("listener")) sim.setListener(params.listenerDist, params.listenerAzimuth, params.listenerElevation);
439+ };
440+ buildControls();
441+ buildViewToolbar();
442+ buildMovieToolbar();
443+ buildSpecToolbar();
444+ const q = new URLSearchParams(location.search);
445+ const set = form.applyUrl(q);
446+ if (set.includes("duration")) {
447+ duration = params.duration;
448+ durInput.value = String(duration);
449+ }
450+ exposeHooks();
451+ if (!("gpu" in navigator)) {
452+ showNoGpu("WebGPU is not available in this browser. Use a recent Chrome or Edge, or Firefox or Safari with WebGPU enabled.");
453+ return;
454+ }
455+ const adapter = (await navigator.gpu.requestAdapter({ powerPreference: "high-performance" })) ?? (await navigator.gpu.requestAdapter());
456+ if (!adapter) {
457+ showNoGpu("No WebGPU adapter was found.");
458+ return;
459+ }
460+ device = await adapter.requestDevice();
461+ device.lost.then((info) => showNoGpu(`WebGPU device lost: ${info.message}`));
462+ device.addEventListener("uncapturederror", (ev) => console.error("WebGPU error:", (ev as GPUUncapturedErrorEvent).error.message));
463+ const { grid } = form.gridInfo();
464+ $("info").innerHTML = `<div class="row"><span>Grid</span><span>${grid.nx} × ${grid.ny} nodes</span></div><p class="sub">Press Start to run.</p>`;
465+ requestAnimationFrame(frame);
466+ if (q.get("selftest") === "1") await selfTest();
467+ else if (q.get("autostart") === "1") void start();
468+}
469+
470+/** Hooks for automated tests, on window.steelsim. */
471+function exposeHooks(): void {
472+ (window as unknown as { steelsim: unknown }).steelsim = {
473+ status: () => (sim ? { ...sim.status(), dt: sim.dt, stepsPerSample: sim.stepsPerSample, nx: sim.grid.nx, ny: sim.grid.ny, stepsPerSec } : null),
474+ audio: () => (sim ? Array.from(sim.recordedAudio()) : null),
475+ frames: () => sim?.frames.length ?? 0,
476+ selftest: () => selfTestResult,
477+ noGpu: () => { const e = $("nogpu"); return e.hidden ? null : e.textContent; },
478+ };
479+}
480+
481+function showNoGpu(msg: string): void {
482+ const el = $("nogpu");
483+ el.hidden = false;
484+ el.textContent = msg;
485+ btnStart.disabled = true;
486+ selfTestResult = { done: true, passed: false, error: msg };
487+}
488+
489+void init();
src/params.tsadded+338−0View file
@@ -0,0 +1,338 @@
1+// Physical and numerical parameters of the simulation, and the derived
2+// quantities (grid, material constants, geometry of the excitation).
3+// All values are SI.
4+
5+export type PlateShape = "rectangle" | "ellipse";
6+export type BowShape = "none" | "width" | "height" | "dish";
7+export type Excitation = "push" | "strike";
8+
9+export interface PhysParams {
10+ // Sheet
11+ shape: PlateShape; // outline of the sheet within its bounding box
12+ Lx: number; // width along x [m] (bounding box for an ellipse)
13+ Ly: number; // height along y [m]
14+ h: number; // thickness [m]
15+ E: number; // Young's modulus [Pa]
16+ nu: number; // Poisson ratio
17+ rho: number; // density [kg/m^3]
18+ holeInset: number; // distance of the thread hole from the edges near the top-left corner [m]
19+ // Stress-free initial shape
20+ bowShape: BowShape;
21+ bowAmplitude: number; // sag of the initial shape at the center [m]
22+ // Losses
23+ sigma0: number; // frequency-independent loss of all motion [1/s]
24+ sigma1: number; // loss proportional to the squared wavenumber, transverse motion [m^2/s]
25+ sigmaIn: number; // additional frequency-independent loss of the in-plane motion [1/s]
26+ // Excitation
27+ excitation: Excitation;
28+ pushAmount: number; // corner displacement [m]; negative pushes against the bow
29+ fingerRadius: number; // radius of the pushed patch [m]
30+ rampTime: number; // duration of the push ramp [s]
31+ holdMax: number; // maximum settle time before the release [s]
32+ releaseTime: number; // time over which the finger force falls to zero [s]
33+ settleSigma0: number; // artificial damping during push and settle [1/s]
34+ settleKE: number; // non-rigid kinetic energy below which the sheet counts as settled [J]
35+ strikeFracX: number; // strike point as a fraction of the width from the left edge
36+ strikeFracY: number; // strike point as a fraction of the height from the bottom edge
37+ strikeVelocity: number; // stick tip speed at impact [m/s]
38+ stickMass: number; // effective mass of the stick [kg]
39+ contactStiffness: number; // Hertz contact stiffness K [N/m^alpha]
40+ contactExponent: number; // Hertz exponent alpha (1.5 for a rounded tip)
41+ contactRadius: number; // radius of the contact patch [m]
42+ // Listener
43+ listenerDist: number; // [m]
44+ listenerAzimuth: number; // [deg], rotation about the sheet's vertical axis (0 = on the normal)
45+ listenerElevation: number; // [deg]
46+ // Numerics and output
47+ spacing: number; // requested grid spacing [m]
48+ fs: number; // audio sample rate [Hz]
49+ fps: number; // movie frame rate, also the batch rate [Hz]
50+ duration: number; // recorded duration after the release or the impact [s]
51+}
52+
53+export const defaultParams: PhysParams = {
54+ shape: "rectangle",
55+ Lx: 12 * 0.0254,
56+ Ly: 18 * 0.0254,
57+ h: 0.79e-3,
58+ E: 193e9,
59+ nu: 0.29,
60+ rho: 7900,
61+ holeInset: 0.012,
62+ bowShape: "width",
63+ bowAmplitude: 3e-3,
64+ sigma0: 1.15,
65+ sigma1: 5e-4,
66+ sigmaIn: 20,
67+ excitation: "push",
68+ pushAmount: 0.01,
69+ fingerRadius: 0.01,
70+ rampTime: 0.1,
71+ holdMax: 0.5,
72+ releaseTime: 0.02,
73+ settleSigma0: 300,
74+ settleKE: 1e-8,
75+ strikeFracX: 0.7,
76+ strikeFracY: 0.35,
77+ strikeVelocity: 3,
78+ stickMass: 0.03,
79+ contactStiffness: 5e8,
80+ contactExponent: 1.5,
81+ contactRadius: 5e-3,
82+ listenerDist: 1.0,
83+ listenerAzimuth: 0,
84+ listenerElevation: 0,
85+ spacing: 0.008,
86+ fs: 48000,
87+ fps: 120,
88+ duration: 4,
89+};
90+
91+export const SPEED_OF_SOUND = 343; // [m/s]
92+export const RHO_AIR = 1.2; // [kg/m^3]
93+export const HIST_LEN = 256; // ring buffer length of the retarded-time listener sum [samples]
94+export const RB_INTERVAL = 32; // steps between rigid-body projections while ringing
95+export const RB_INTERVAL_HELD = 4; // while the finger constraint is active
96+
97+// ---------- Grid ----------
98+
99+export interface Grid {
100+ nx: number;
101+ ny: number;
102+ hx: number;
103+ hy: number;
104+ n: number;
105+ /** 1 for nodes that belong to the sheet, 0 otherwise; absent for a full rectangle. */
106+ mask?: Uint8Array;
107+}
108+
109+export function makeGrid(p: PhysParams): Grid {
110+ const nx = Math.max(6, Math.round(p.Lx / p.spacing) + 1);
111+ const ny = Math.max(6, Math.round(p.Ly / p.spacing) + 1);
112+ const g: Grid = { nx, ny, hx: p.Lx / (nx - 1), hy: p.Ly / (ny - 1), n: nx * ny };
113+ return p.shape === "ellipse" ? withEllipseMask(g) : g;
114+}
115+
116+/**
117+ * Restricts a grid to the ellipse inscribed in its bounding box. The edge is a
118+ * staircase of whole cells: a cell belongs to the sheet when its center lies
119+ * inside the ellipse, and a node when it is a corner of such a cell. Deciding
120+ * by cell centers rather than by nodes keeps the covered area unbiased; a
121+ * node-based mask loses about half a cell all around the edge, which raises
122+ * every frequency by a few percent. A node with no cell of its own would carry
123+ * no mass, and the trailing loop removes any such node.
124+ */
125+export function withEllipseMask(g: Grid): Grid {
126+ const mask = new Uint8Array(g.n);
127+ const a = ((g.nx - 1) * g.hx) / 2, b = ((g.ny - 1) * g.hy) / 2;
128+ for (let cj = 0; cj < g.ny - 1; cj++) {
129+ for (let ci = 0; ci < g.nx - 1; ci++) {
130+ const X = (nodeX(g, ci) + g.hx / 2) / a, Y = (nodeY(g, cj) + g.hy / 2) / b;
131+ if (X * X + Y * Y <= 1) {
132+ mask[cj * g.nx + ci] = 1; mask[cj * g.nx + ci + 1] = 1;
133+ mask[(cj + 1) * g.nx + ci] = 1; mask[(cj + 1) * g.nx + ci + 1] = 1;
134+ }
135+ }
136+ }
137+ const gm: Grid = { ...g, mask };
138+ for (let changed = true; changed; ) {
139+ changed = false;
140+ for (let j = 0; j < g.ny; j++) {
141+ for (let i = 0; i < g.nx; i++) {
142+ if (mask[j * g.nx + i] && nodeWeight(gm, i, j) === 0) {
143+ mask[j * g.nx + i] = 0;
144+ changed = true;
145+ }
146+ }
147+ }
148+ }
149+ return gm;
150+}
151+
152+/** Node coordinates relative to the center of the bounding box. */
153+export function nodeX(g: Grid, i: number): number {
154+ return (i - (g.nx - 1) / 2) * g.hx;
155+}
156+export function nodeY(g: Grid, j: number): number {
157+ return (j - (g.ny - 1) / 2) * g.hy;
158+}
159+
160+/** True when node (i, j) exists and belongs to the sheet. */
161+export function insideNode(g: Grid, i: number, j: number): boolean {
162+ if (i < 0 || i >= g.nx || j < 0 || j >= g.ny) return false;
163+ return g.mask ? g.mask[j * g.nx + i] !== 0 : true;
164+}
165+
166+/** True when cell (ci, cj), with nodes (ci..ci+1, cj..cj+1), lies within the sheet. */
167+export function insideCell(g: Grid, ci: number, cj: number): boolean {
168+ return insideNode(g, ci, cj) && insideNode(g, ci + 1, cj) && insideNode(g, ci + 1, cj + 1) && insideNode(g, ci, cj + 1);
169+}
170+
171+/**
172+ * Quadrature and lumped-mass weight of a node: the fraction of its four
173+ * surrounding cells that belong to the sheet. On a full rectangle this is the
174+ * trapezoidal rule (1 inside, 1/2 on an edge, 1/4 at a corner). Zero outside.
175+ */
176+export function nodeWeight(g: Grid, i: number, j: number): number {
177+ if (!insideNode(g, i, j)) return 0;
178+ let c = 0;
179+ if (insideCell(g, i, j)) c++;
180+ if (insideCell(g, i - 1, j)) c++;
181+ if (insideCell(g, i, j - 1)) c++;
182+ if (insideCell(g, i - 1, j - 1)) c++;
183+ return c / 4;
184+}
185+
186+/** Per-node weights as a flat array. */
187+export function nodeWeights(g: Grid): Float32Array<ArrayBuffer> {
188+ const w = new Float32Array(g.n);
189+ for (let j = 0; j < g.ny; j++) for (let i = 0; i < g.nx; i++) w[j * g.nx + i] = nodeWeight(g, i, j);
190+ return w;
191+}
192+
193+/** Number of nodes in the sheet. */
194+export function countInside(g: Grid): number {
195+ if (!g.mask) return g.n;
196+ let c = 0;
197+ for (let k = 0; k < g.n; k++) c += g.mask[k];
198+ return c;
199+}
200+
201+// ---------- Material ----------
202+
203+export interface Material {
204+ D: number; // bending stiffness E h^3 / (12 (1 - nu^2)) [N m]
205+ C: number; // membrane stiffness E h / (1 - nu^2) [N/m]
206+ G: number; // membrane shear stiffness E h / (2 (1 + nu)) [N/m]
207+ nu: number;
208+ rhoh: number; // mass per unit area [kg/m^2]
209+}
210+
211+export function makeMaterial(p: PhysParams): Material {
212+ return {
213+ D: (p.E * p.h ** 3) / (12 * (1 - p.nu * p.nu)),
214+ C: (p.E * p.h) / (1 - p.nu * p.nu),
215+ G: (p.E * p.h) / (2 * (1 + p.nu)),
216+ nu: p.nu,
217+ rhoh: p.rho * p.h,
218+ };
219+}
220+
221+/**
222+ * 60 dB decay time [s] of a bending mode of angular frequency omega under the
223+ * loss model sigma0 + sigma1 k^2, using the flat-plate dispersion k^2 = omega / kappa.
224+ */
225+export function t60(m: Material, p: PhysParams, omega: number): number {
226+ const kappa = Math.sqrt(m.D / m.rhoh);
227+ const sigma = p.sigma0 + p.sigma1 * (omega / kappa);
228+ return (3 * Math.LN10) / sigma;
229+}
230+
231+/**
232+ * Frequency [Hz] at which the three-point second difference lowers a bending
233+ * frequency by the relative amount `err` on a grid of spacing hmax. A plane
234+ * bending wave along a grid axis with wavenumber k has the discrete frequency
235+ * ratio (2 - 2 cos kh) / (kh)^2.
236+ */
237+export function dispersionFrequency(m: Material, g: Grid, err: number): number {
238+ const h = Math.max(g.hx, g.hy);
239+ let lo = 0, hi = Math.PI;
240+ for (let it = 0; it < 60; it++) {
241+ const x = (lo + hi) / 2;
242+ if ((2 - 2 * Math.cos(x)) / (x * x) > 1 - err) lo = x; else hi = x;
243+ }
244+ const k = lo / h;
245+ const kappa = Math.sqrt(m.D / m.rhoh);
246+ return (kappa * k * k) / (2 * Math.PI);
247+}
248+
249+// ---------- Geometry of the excitation ----------
250+
251+/** Strike point in sheet coordinates. */
252+export function strikePoint(p: PhysParams): [number, number] {
253+ return [(p.strikeFracX - 0.5) * p.Lx, (p.strikeFracY - 0.5) * p.Ly];
254+}
255+
256+/**
257+ * Center of the thread hole. Rectangle: inset from the two edges of the
258+ * top-left corner. Ellipse: inset along the radius from the rim point in the
259+ * top-left direction (parametric angle 135 degrees).
260+ */
261+export function holePoint(p: PhysParams): [number, number] {
262+ if (p.shape === "ellipse") {
263+ const c = Math.SQRT1_2;
264+ return [-(p.Lx / 2 - p.holeInset) * c, (p.Ly / 2 - p.holeInset) * c];
265+ }
266+ return [-p.Lx / 2 + p.holeInset, p.Ly / 2 - p.holeInset];
267+}
268+
269+/** Center of the pushed patch: the bottom-right corner, or the rim point opposite the hole. */
270+export function pushPoint(p: PhysParams): [number, number] {
271+ if (p.shape === "ellipse") {
272+ const c = Math.SQRT1_2;
273+ return [(p.Lx / 2) * c, -(p.Ly / 2) * c];
274+ }
275+ return [p.Lx / 2, -p.Ly / 2];
276+}
277+
278+/**
279+ * Stress-free initial shape w0 at every node. All three shapes are parabolic
280+ * with the given sag at the center of the bounding box: "width" and "height"
281+ * are cylindrical with zero sag along the two opposite edges, and "dish" has
282+ * the same curvature in both directions, so its sag is measured relative to
283+ * the corners (the edge midpoints sit at half the sag).
284+ */
285+export function initialShape(p: PhysParams, g: Grid): Float32Array<ArrayBuffer> {
286+ const w0 = new Float32Array(g.n);
287+ if (p.bowShape === "none" || p.bowAmplitude === 0) return w0;
288+ for (let j = 0; j < g.ny; j++) {
289+ for (let i = 0; i < g.nx; i++) {
290+ const X = (2 * nodeX(g, i)) / p.Lx; // -1..1
291+ const Y = (2 * nodeY(g, j)) / p.Ly;
292+ let f = 0;
293+ if (p.bowShape === "width") f = 1 - X * X;
294+ else if (p.bowShape === "height") f = 1 - Y * Y;
295+ else f = 1 - 0.5 * (X * X + Y * Y);
296+ w0[j * g.nx + i] = p.bowAmplitude * f;
297+ }
298+ }
299+ return w0;
300+}
301+
302+/**
303+ * Nodes of the drumstick contact patch: the sheet node nearest the strike
304+ * point, plus every sheet node within the contact radius. The contact force
305+ * is spread over these nodes in proportion to their mass, and the contact
306+ * coordinate is their mass-weighted mean displacement, so the two are
307+ * conjugate and the contact has a potential energy.
308+ */
309+export function contactPatch(p: PhysParams, g: Grid): { nodes: number[]; nearest: number; mass: number; box: [number, number, number, number] } {
310+ const m = makeMaterial(p);
311+ const [sx, sy] = strikePoint(p);
312+ let nearest = -1, best = Infinity;
313+ for (let j = 0; j < g.ny; j++) {
314+ for (let i = 0; i < g.nx; i++) {
315+ if (nodeWeight(g, i, j) === 0) continue;
316+ const d2 = (nodeX(g, i) - sx) ** 2 + (nodeY(g, j) - sy) ** 2;
317+ if (d2 < best) { best = d2; nearest = j * g.nx + i; }
318+ }
319+ }
320+ const nodes: number[] = [];
321+ let mass = 0;
322+ const box: [number, number, number, number] = [g.nx, -1, g.ny, -1]; // i0, i1, j0, j1 (inclusive)
323+ for (let j = 0; j < g.ny; j++) {
324+ for (let i = 0; i < g.nx; i++) {
325+ const k = j * g.nx + i;
326+ const wgt = nodeWeight(g, i, j);
327+ if (wgt === 0) continue;
328+ const d2 = (nodeX(g, i) - sx) ** 2 + (nodeY(g, j) - sy) ** 2;
329+ if (d2 <= p.contactRadius ** 2 || k === nearest) {
330+ nodes.push(k);
331+ mass += m.rhoh * g.hx * g.hy * wgt;
332+ box[0] = Math.min(box[0], i); box[1] = Math.max(box[1], i);
333+ box[2] = Math.min(box[2], j); box[3] = Math.max(box[3], j);
334+ }
335+ }
336+ }
337+ return { nodes, nearest, mass, box };
338+}
src/render/mat4.tsadded+50−0View file
@@ -0,0 +1,50 @@
1+// Minimal column-major 4x4 matrix helpers (WebGPU layout).
2+
3+export type Mat4 = Float32Array;
4+
5+export function perspective(fovY: number, aspect: number, near: number, far: number): Mat4 {
6+ const f = 1 / Math.tan(fovY / 2);
7+ const m = new Float32Array(16);
8+ m[0] = f / aspect;
9+ m[5] = f;
10+ m[10] = far / (near - far);
11+ m[11] = -1;
12+ m[14] = (near * far) / (near - far);
13+ return m;
14+}
15+
16+export function lookAt(eye: number[], target: number[], up: number[]): Mat4 {
17+ const z = normalize([eye[0] - target[0], eye[1] - target[1], eye[2] - target[2]]);
18+ const x = normalize(cross(up, z));
19+ const y = cross(z, x);
20+ const m = new Float32Array(16);
21+ m[0] = x[0]; m[1] = y[0]; m[2] = z[0];
22+ m[4] = x[1]; m[5] = y[1]; m[6] = z[1];
23+ m[8] = x[2]; m[9] = y[2]; m[10] = z[2];
24+ m[12] = -dot(x, eye); m[13] = -dot(y, eye); m[14] = -dot(z, eye);
25+ m[15] = 1;
26+ return m;
27+}
28+
29+export function multiply(a: Mat4, b: Mat4): Mat4 {
30+ const out = new Float32Array(16);
31+ for (let c = 0; c < 4; c++) {
32+ for (let r = 0; r < 4; r++) {
33+ let s = 0;
34+ for (let k = 0; k < 4; k++) s += a[k * 4 + r] * b[c * 4 + k];
35+ out[c * 4 + r] = s;
36+ }
37+ }
38+ return out;
39+}
40+
41+function cross(a: number[], b: number[]): number[] {
42+ return [a[1] * b[2] - a[2] * b[1], a[2] * b[0] - a[0] * b[2], a[0] * b[1] - a[1] * b[0]];
43+}
44+function dot(a: number[], b: number[]): number {
45+ return a[0] * b[0] + a[1] * b[1] + a[2] * b[2];
46+}
47+function normalize(a: number[]): number[] {
48+ const l = Math.hypot(a[0], a[1], a[2]) || 1;
49+ return [a[0] / l, a[1] / l, a[2] / l];
50+}
src/render/renderer.tsadded+187−0View file
@@ -0,0 +1,187 @@
1+// WebGPU rendering of the displacement field: a 2D color map or a lit 3D height field.
2+
3+import { type Grid, type PhysParams, holePoint, insideCell, pushPoint, strikePoint } from "../params.ts";
4+import { lookAt, multiply, perspective } from "./mat4.ts";
5+import shaderCode from "./shaders.wgsl?raw";
6+
7+export type ViewMode = "2d" | "3d";
8+
9+export class Renderer {
10+ mode: ViewMode = "2d";
11+ /** Vertical exaggeration of the 3D view. */
12+ exaggeration = 4;
13+ /** Color range |w| [m]; 0 means automatic (from the current maximum). */
14+ colorRange = 0;
15+ /** Add the stress-free initial shape to the 3D height field. */
16+ showShape = false;
17+ camera = { yaw: 0.6, pitch: 0.5, dist: 1.0 };
18+
19+ private device: GPUDevice;
20+ private canvas: HTMLCanvasElement;
21+ private ctx: GPUCanvasContext;
22+ private format: GPUTextureFormat;
23+ private uniform: GPUBuffer;
24+ private pipe2d: GPURenderPipeline;
25+ private pipe3d: GPURenderPipeline;
26+ private bgLayout: GPUBindGroupLayout;
27+ private index: GPUBuffer;
28+ private indexCount: number;
29+ private depth: GPUTexture | null = null;
30+ private grid: Grid;
31+ private params: PhysParams;
32+ private bindGroups = new Map<GPUBuffer, GPUBindGroup>();
33+
34+ private w0: GPUBuffer;
35+ private nw: GPUBuffer;
36+
37+ constructor(device: GPUDevice, canvas: HTMLCanvasElement, grid: Grid, params: PhysParams, w0: GPUBuffer, nw: GPUBuffer) {
38+ this.device = device;
39+ this.w0 = w0;
40+ this.nw = nw;
41+ this.canvas = canvas;
42+ this.grid = grid;
43+ this.params = params;
44+ this.ctx = canvas.getContext("webgpu") as GPUCanvasContext;
45+ this.format = navigator.gpu.getPreferredCanvasFormat();
46+ this.ctx.configure({ device, format: this.format, alphaMode: "opaque" });
47+ this.uniform = device.createBuffer({ size: 144, usage: GPUBufferUsage.UNIFORM | GPUBufferUsage.COPY_DST });
48+ const module = device.createShaderModule({ code: shaderCode });
49+ this.bgLayout = device.createBindGroupLayout({
50+ entries: [
51+ { binding: 0, visibility: GPUShaderStage.VERTEX | GPUShaderStage.FRAGMENT, buffer: { type: "uniform" } },
52+ { binding: 1, visibility: GPUShaderStage.VERTEX | GPUShaderStage.FRAGMENT, buffer: { type: "read-only-storage" } },
53+ { binding: 2, visibility: GPUShaderStage.VERTEX | GPUShaderStage.FRAGMENT, buffer: { type: "read-only-storage" } },
54+ { binding: 3, visibility: GPUShaderStage.VERTEX | GPUShaderStage.FRAGMENT, buffer: { type: "read-only-storage" } },
55+ ],
56+ });
57+ const layout = device.createPipelineLayout({ bindGroupLayouts: [this.bgLayout] });
58+ this.pipe2d = device.createRenderPipeline({
59+ layout,
60+ vertex: { module, entryPoint: "vs2d" },
61+ fragment: { module, entryPoint: "fs2d", targets: [{ format: this.format }] },
62+ primitive: { topology: "triangle-list" },
63+ });
64+ this.pipe3d = device.createRenderPipeline({
65+ layout,
66+ vertex: { module, entryPoint: "vs3d" },
67+ fragment: { module, entryPoint: "fs3d", targets: [{ format: this.format }] },
68+ primitive: { topology: "triangle-list", cullMode: "none" },
69+ depthStencil: { format: "depth24plus", depthWriteEnabled: true, depthCompare: "less" },
70+ });
71+ // Triangle indices of the height field (cells of the sheet only).
72+ const { nx, ny } = grid;
73+ const idx = new Uint32Array((nx - 1) * (ny - 1) * 6);
74+ let k = 0;
75+ for (let j = 0; j < ny - 1; j++) {
76+ for (let i = 0; i < nx - 1; i++) {
77+ if (!insideCell(grid, i, j)) continue;
78+ const a = j * nx + i, b = a + 1, c = a + nx, d = c + 1;
79+ idx[k++] = a; idx[k++] = b; idx[k++] = d;
80+ idx[k++] = a; idx[k++] = d; idx[k++] = c;
81+ }
82+ }
83+ this.indexCount = k;
84+ this.index = device.createBuffer({ size: Math.max(16, idx.byteLength), usage: GPUBufferUsage.INDEX | GPUBufferUsage.COPY_DST });
85+ device.queue.writeBuffer(this.index, 0, idx);
86+ this.installControls();
87+ }
88+
89+ private installControls(): void {
90+ let dragging = false, lastX = 0, lastY = 0;
91+ this.canvas.addEventListener("pointerdown", (e) => {
92+ dragging = true; lastX = e.clientX; lastY = e.clientY;
93+ this.canvas.setPointerCapture(e.pointerId);
94+ });
95+ this.canvas.addEventListener("pointerup", () => (dragging = false));
96+ this.canvas.addEventListener("pointermove", (e) => {
97+ if (!dragging || this.mode !== "3d") return;
98+ this.camera.yaw += (e.clientX - lastX) * 0.01;
99+ this.camera.pitch = Math.max(-1.5, Math.min(1.5, this.camera.pitch + (e.clientY - lastY) * 0.01));
100+ lastX = e.clientX; lastY = e.clientY;
101+ });
102+ this.canvas.addEventListener("wheel", (e) => {
103+ if (this.mode !== "3d") return;
104+ e.preventDefault();
105+ this.camera.dist *= Math.exp(e.deltaY * 0.001);
106+ this.camera.dist = Math.max(0.2, Math.min(5, this.camera.dist));
107+ }, { passive: false });
108+ }
109+
110+ private resize(): void {
111+ const dpr = window.devicePixelRatio || 1;
112+ const w = Math.max(1, Math.floor(this.canvas.clientWidth * dpr));
113+ const h = Math.max(1, Math.floor(this.canvas.clientHeight * dpr));
114+ if (this.canvas.width !== w || this.canvas.height !== h || !this.depth) {
115+ this.canvas.width = w;
116+ this.canvas.height = h;
117+ this.depth?.destroy();
118+ this.depth = this.device.createTexture({ size: [w, h], format: "depth24plus", usage: GPUTextureUsage.RENDER_ATTACHMENT });
119+ }
120+ }
121+
122+ /** Renders the w field held in `wBuffer` (f32 per node). `wmax` drives the automatic color range. */
123+ render(encoder: GPUCommandEncoder, wBuffer: GPUBuffer, wmax: number): void {
124+ this.resize();
125+ const aspect = this.canvas.width / this.canvas.height;
126+ const range = this.colorRange > 0 ? this.colorRange : Math.max(wmax, 1e-6);
127+ const u = new Float32Array(36);
128+ const ui = new Uint32Array(u.buffer);
129+ const eye = [
130+ this.camera.dist * Math.cos(this.camera.pitch) * Math.sin(this.camera.yaw),
131+ this.camera.dist * Math.sin(this.camera.pitch),
132+ this.camera.dist * Math.cos(this.camera.pitch) * Math.cos(this.camera.yaw),
133+ ];
134+ const mvp = multiply(perspective(0.7, aspect, 0.02, 20), lookAt(eye, [0, 0, 0], [0, 1, 0]));
135+ u.set(mvp, 0);
136+ ui[16] = this.grid.nx; ui[17] = this.grid.ny;
137+ u[18] = this.grid.hx; u[19] = this.grid.hy;
138+ u[20] = this.exaggeration; u[21] = range; u[22] = aspect;
139+ ui[23] = this.mode === "3d" ? 1 : 0;
140+ const p = this.params;
141+ [u[24], u[25]] = holePoint(p);
142+ if (p.excitation === "strike") {
143+ const [sx, sy] = strikePoint(p);
144+ u[26] = sx; u[27] = sy; u[28] = p.contactRadius;
145+ } else {
146+ [u[26], u[27]] = pushPoint(p);
147+ u[28] = p.fingerRadius;
148+ }
149+ u[29] = 0.4; u[30] = 0.5; u[31] = 1.0;
150+ ui[32] = this.showShape ? 1 : 0;
151+ ui[33] = p.shape === "ellipse" ? 1 : 0;
152+ this.device.queue.writeBuffer(this.uniform, 0, u);
153+
154+ let bg = this.bindGroups.get(wBuffer);
155+ if (!bg) {
156+ bg = this.device.createBindGroup({
157+ layout: this.bgLayout,
158+ entries: [
159+ { binding: 0, resource: { buffer: this.uniform } },
160+ { binding: 1, resource: { buffer: wBuffer } },
161+ { binding: 2, resource: { buffer: this.w0 } },
162+ { binding: 3, resource: { buffer: this.nw } },
163+ ],
164+ });
165+ this.bindGroups.set(wBuffer, bg);
166+ }
167+ const view = this.ctx.getCurrentTexture().createView();
168+ const clear = { r: 0.11, g: 0.12, b: 0.14, a: 1 };
169+ if (this.mode === "2d") {
170+ const pass = encoder.beginRenderPass({ colorAttachments: [{ view, clearValue: clear, loadOp: "clear", storeOp: "store" }] });
171+ pass.setPipeline(this.pipe2d);
172+ pass.setBindGroup(0, bg);
173+ pass.draw(3);
174+ pass.end();
175+ } else {
176+ const pass = encoder.beginRenderPass({
177+ colorAttachments: [{ view, clearValue: clear, loadOp: "clear", storeOp: "store" }],
178+ depthStencilAttachment: { view: this.depth!.createView(), depthClearValue: 1, depthLoadOp: "clear", depthStoreOp: "store" },
179+ });
180+ pass.setPipeline(this.pipe3d);
181+ pass.setBindGroup(0, bg);
182+ pass.setIndexBuffer(this.index, "uint32");
183+ pass.drawIndexed(this.indexCount);
184+ pass.end();
185+ }
186+ }
187+}
src/render/shaders.wgsladded+145−0View file
@@ -0,0 +1,145 @@
1+struct RU {
2+ mvp: mat4x4f,
3+ nx: u32,
4+ ny: u32,
5+ hx: f32,
6+ hy: f32,
7+ wScale: f32, // vertical exaggeration in the 3D view
8+ colorRange: f32, // |w| mapped to full color [m]
9+ aspect: f32, // canvas width / height
10+ mode: u32, // 0: 2D, 1: 3D
11+ holeX: f32,
12+ holeY: f32,
13+ pushX: f32,
14+ pushY: f32,
15+ markerR: f32,
16+ lightX: f32,
17+ lightY: f32,
18+ lightZ: f32,
19+ showShape: u32, // 3D: add the stress-free initial shape to the height field
20+ shape: u32, // 0: rectangle, 1: ellipse inscribed in the bounding box
21+ pad1: u32,
22+ pad2: u32,
23+};
24+
25+@group(0) @binding(0) var<uniform> R: RU;
26+@group(0) @binding(1) var<storage, read> w: array<f32>;
27+@group(0) @binding(2) var<storage, read> w0: array<f32>;
28+@group(0) @binding(3) var<storage, read> nw: array<f32>; // node weights, 0 outside the sheet
29+
30+fn colormap(t: f32) -> vec3f {
31+ let s = clamp(t, -1.0, 1.0);
32+ let mid = vec3f(0.93, 0.93, 0.93);
33+ let blue = vec3f(0.10, 0.35, 0.85);
34+ let red = vec3f(0.85, 0.20, 0.12);
35+ if (s < 0.0) { return mix(mid, blue, -s); }
36+ return mix(mid, red, s);
37+}
38+
39+fn wAt(i: i32, j: i32) -> f32 {
40+ let ii = clamp(i, 0, i32(R.nx) - 1);
41+ let jj = clamp(j, 0, i32(R.ny) - 1);
42+ return w[u32(jj) * R.nx + u32(ii)];
43+}
44+fn nwAt(i: i32, j: i32) -> f32 {
45+ if (i < 0 || i >= i32(R.nx) || j < 0 || j >= i32(R.ny)) { return 0.0; }
46+ return nw[u32(j) * R.nx + u32(i)];
47+}
48+// Height used by the 3D view: w, plus the initial shape when requested.
49+fn zAt(i: i32, j: i32) -> f32 {
50+ let ii = clamp(i, 0, i32(R.nx) - 1);
51+ let jj = clamp(j, 0, i32(R.ny) - 1);
52+ var z = w[u32(jj) * R.nx + u32(ii)];
53+ if (R.showShape != 0u) { z += w0[u32(jj) * R.nx + u32(ii)]; }
54+ return z;
55+}
56+// Slope along one axis, one-sided where a neighbor lies outside the sheet.
57+fn slope(i: i32, j: i32, di: i32, dj: i32, h: f32) -> f32 {
58+ let hasP = nwAt(i + di, j + dj) > 0.0;
59+ let hasM = nwAt(i - di, j - dj) > 0.0;
60+ if (hasP && hasM) { return (zAt(i + di, j + dj) - zAt(i - di, j - dj)) / (2.0 * h); }
61+ if (hasP) { return (zAt(i + di, j + dj) - zAt(i, j)) / h; }
62+ if (hasM) { return (zAt(i, j) - zAt(i - di, j - dj)) / h; }
63+ return 0.0;
64+}
65+
66+// ---------- 2D view: full-screen triangle, bilinear lookup in the fragment shader ----------
67+
68+struct V2 { @builtin(position) pos: vec4f, @location(0) uv: vec2f };
69+
70+@vertex
71+fn vs2d(@builtin(vertex_index) vi: u32) -> V2 {
72+ var p = array<vec2f, 3>(vec2f(-1.0, -1.0), vec2f(3.0, -1.0), vec2f(-1.0, 3.0));
73+ var o: V2;
74+ o.pos = vec4f(p[vi], 0.0, 1.0);
75+ o.uv = p[vi];
76+ return o;
77+}
78+
79+@fragment
80+fn fs2d(in: V2) -> @location(0) vec4f {
81+ let Lx = f32(R.nx - 1u) * R.hx;
82+ let Ly = f32(R.ny - 1u) * R.hy;
83+ // Fit the plate into the canvas with a margin, preserving its aspect ratio.
84+ let margin = 0.92;
85+ var sx = margin; var sy = margin;
86+ let plateAspect = Lx / Ly;
87+ if (plateAspect > R.aspect) { sy = margin * R.aspect / plateAspect; } else { sx = margin * plateAspect / R.aspect; }
88+ let px = in.uv.x / sx * 0.5 * Lx;
89+ let py = in.uv.y / sy * 0.5 * Ly;
90+ let bg = vec4f(0.11, 0.12, 0.14, 1.0);
91+ if (abs(px) > 0.5 * Lx || abs(py) > 0.5 * Ly) { return bg; }
92+ if (R.shape == 1u) {
93+ let ex = px / (0.5 * Lx); let ey = py / (0.5 * Ly);
94+ if (ex * ex + ey * ey > 1.0) { return bg; }
95+ }
96+ let fi = (px + 0.5 * Lx) / R.hx;
97+ let fj = (py + 0.5 * Ly) / R.hy;
98+ let i0 = i32(floor(fi)); let j0 = i32(floor(fj));
99+ let a = fi - f32(i0); let b = fj - f32(j0);
100+ // Bilinear interpolation over the nodes of the sheet only (outside nodes carry no value).
101+ let c = vec4f((1.0 - a) * (1.0 - b), a * (1.0 - b), (1.0 - a) * b, a * b);
102+ let m = vec4f(nwAt(i0, j0), nwAt(i0 + 1, j0), nwAt(i0, j0 + 1), nwAt(i0 + 1, j0 + 1));
103+ let wgt = c * step(vec4f(1e-6), m);
104+ let den = wgt.x + wgt.y + wgt.z + wgt.w;
105+ if (den <= 0.0) { return bg; }
106+ let val = dot(wgt, vec4f(wAt(i0, j0), wAt(i0 + 1, j0), wAt(i0, j0 + 1), wAt(i0 + 1, j0 + 1))) / den;
107+ var col = colormap(val / R.colorRange);
108+ // Markers: hole (ring) and pushed patch (dashed-looking ring).
109+ let dh = distance(vec2f(px, py), vec2f(R.holeX, R.holeY));
110+ if (abs(dh - R.markerR) < 0.35 * R.markerR) { col = vec3f(0.1, 0.1, 0.1); }
111+ if (dh < 0.4 * R.markerR) { col = vec3f(0.11, 0.12, 0.14); }
112+ let dp = distance(vec2f(px, py), vec2f(R.pushX, R.pushY));
113+ if (abs(dp - R.markerR * 1.6) < 0.25 * R.markerR) { col = vec3f(0.95, 0.75, 0.1); }
114+ return vec4f(col, 1.0);
115+}
116+
117+// ---------- 3D view: height field ----------
118+
119+struct V3 { @builtin(position) pos: vec4f, @location(0) val: f32, @location(1) normal: vec3f };
120+
121+@vertex
122+fn vs3d(@builtin(vertex_index) vi: u32) -> V3 {
123+ let i = i32(vi % R.nx); let j = i32(vi / R.nx);
124+ let x = (f32(i) - 0.5 * f32(R.nx - 1u)) * R.hx;
125+ let y = (f32(j) - 0.5 * f32(R.ny - 1u)) * R.hy;
126+ let val = wAt(i, j);
127+ let z = zAt(i, j) * R.wScale;
128+ let dzdx = slope(i, j, 1, 0, R.hx) * R.wScale;
129+ let dzdy = slope(i, j, 0, 1, R.hy) * R.wScale;
130+ var o: V3;
131+ o.pos = R.mvp * vec4f(x, y, z, 1.0);
132+ o.val = val;
133+ o.normal = normalize(vec3f(-dzdx, -dzdy, 1.0));
134+ return o;
135+}
136+
137+@fragment
138+fn fs3d(in: V3) -> @location(0) vec4f {
139+ let l = normalize(vec3f(R.lightX, R.lightY, R.lightZ));
140+ let n = normalize(in.normal);
141+ let diff = abs(dot(n, l));
142+ let base = colormap(in.val / R.colorRange);
143+ let col = base * (0.35 + 0.65 * diff) + vec3f(0.15) * pow(diff, 24.0);
144+ return vec4f(col, 1.0);
145+}
src/selftest.tsadded+111−0View file
@@ -0,0 +1,111 @@
1+// GPU-versus-CPU self-test: one time step of the WGSL kernel on a smooth,
2+// strongly nonlinear random state is compared with the float64 reference.
3+// Runs in the browser (?selftest=1) and is driven by tests/gpu.test.ts.
4+
5+import { defaultParams, initialShape, makeMaterial, nodeWeight, type PhysParams } from "./params.ts";
6+import { accelerations, addDamping, fields } from "./sim/cpuRef.ts";
7+import { Simulation } from "./sim/simulation.ts";
8+
9+export interface SelfTestCase {
10+ name: string;
11+ /** Maximum relative error of the acceleration of each component, over the sheet. */
12+ accError: { u: number; v: number; w: number };
13+ /** Maximum error of the compensated w update relative to the largest increment dt * w_t. */
14+ positionError: number;
15+ nodes: number;
16+ dt: number;
17+}
18+
19+export const ACC_TOLERANCE = 1e-3;
20+export const POSITION_TOLERANCE = 1e-2;
21+
22+const cases: { name: string; overrides: Partial<PhysParams> }[] = [
23+ { name: "rectangle, flat", overrides: { shape: "rectangle", bowShape: "none" } },
24+ { name: "rectangle, bowed across the width", overrides: { shape: "rectangle", bowShape: "width" } },
25+ { name: "ellipse, dish", overrides: { shape: "ellipse", bowShape: "dish" } },
26+];
27+
28+/** A smooth random field: a few sinusoids across the sheet plus a little white noise. */
29+function smoothField(nx: number, ny: number, amplitude: number, noise: number, rnd: () => number): Float64Array {
30+ const f = new Float64Array(nx * ny);
31+ const terms = Array.from({ length: 6 }, () => ({ kx: 1 + Math.floor(4 * rnd()), ky: 1 + Math.floor(4 * rnd()), px: 2 * Math.PI * rnd(), py: 2 * Math.PI * rnd(), a: rnd() - 0.5 }));
32+ for (let j = 0; j < ny; j++) {
33+ for (let i = 0; i < nx; i++) {
34+ let s = 0;
35+ for (const t of terms) s += t.a * Math.cos((t.kx * Math.PI * i) / (nx - 1) + t.px) * Math.cos((t.ky * Math.PI * j) / (ny - 1) + t.py);
36+ f[j * nx + i] = amplitude * s + noise * (rnd() - 0.5);
37+ }
38+ }
39+ return f;
40+}
41+
42+export async function runSelfTest(device: GPUDevice): Promise<SelfTestCase[]> {
43+ const results: SelfTestCase[] = [];
44+ let seed = 2024;
45+ const rnd = () => ((seed = (seed * 1664525 + 1013904223) >>> 0) / 4294967296);
46+ for (const c of cases) {
47+ const params: PhysParams = { ...defaultParams, ...c.overrides, spacing: 0.012 };
48+ const sim = new Simulation(device, params);
49+ const g = sim.grid;
50+ const n = g.n;
51+ const m = makeMaterial(params);
52+ const w0 = Float64Array.from(initialShape(params, g));
53+
54+ // State in float32 (what the GPU sees), then the same values in float64 for the reference.
55+ const pos = new Float32Array(4 * n), vel = new Float32Array(3 * n);
56+ const fw = smoothField(g.nx, g.ny, 2e-3, 5e-5, rnd); // about 2.5 thicknesses: strongly nonlinear
57+ const fu = smoothField(g.nx, g.ny, 2e-5, 1e-6, rnd);
58+ const fv = smoothField(g.nx, g.ny, 2e-5, 1e-6, rnd);
59+ const fdw = smoothField(g.nx, g.ny, 0.5, 1e-2, rnd);
60+ const fdu = smoothField(g.nx, g.ny, 1e-2, 1e-3, rnd);
61+ const fdv = smoothField(g.nx, g.ny, 1e-2, 1e-3, rnd);
62+ for (let k = 0; k < n; k++) {
63+ if (nodeWeight(g, k % g.nx, Math.floor(k / g.nx)) === 0) continue;
64+ pos[k] = fu[k]; pos[n + k] = fv[k]; pos[2 * n + k] = fw[k]; pos[3 * n + k] = 0;
65+ vel[k] = fdu[k]; vel[n + k] = fdv[k]; vel[2 * n + k] = fdw[k];
66+ }
67+ const x = fields(n), v = fields(n);
68+ for (let k = 0; k < n; k++) {
69+ x.u[k] = pos[k]; x.v[k] = pos[n + k]; x.w[k] = pos[2 * n + k];
70+ v.u[k] = vel[k]; v.v[k] = vel[n + k]; v.w[k] = vel[2 * n + k];
71+ }
72+
73+ sim.writeState(pos, vel);
74+ await sim.stepOnce();
75+ const out = await sim.readState();
76+ sim.destroy();
77+
78+ // Reference: a = elastic + losses, v' = v + dt a, w' = w + dt v'.
79+ const acc = fields(n);
80+ accelerations(g, m, x, acc, w0);
81+ addDamping(g, params.sigma0, params.sigma1, params.sigmaIn, v, acc);
82+ const dt = sim.dt;
83+ const err = { u: 0, v: 0, w: 0 };
84+ let scale = { u: 0, v: 0, w: 0 };
85+ let posErr = 0, posScale = 0;
86+ for (let k = 0; k < n; k++) {
87+ if (nodeWeight(g, k % g.nx, Math.floor(k / g.nx)) === 0) continue;
88+ const gu = (out.vel[k] - vel[k]) / dt, gv = (out.vel[n + k] - vel[n + k]) / dt, gw = (out.vel[2 * n + k] - vel[2 * n + k]) / dt;
89+ err.u = Math.max(err.u, Math.abs(gu - acc.u[k]));
90+ err.v = Math.max(err.v, Math.abs(gv - acc.v[k]));
91+ err.w = Math.max(err.w, Math.abs(gw - acc.w[k]));
92+ scale = { u: Math.max(scale.u, Math.abs(acc.u[k])), v: Math.max(scale.v, Math.abs(acc.v[k])), w: Math.max(scale.w, Math.abs(acc.w[k])) };
93+ const vNew = v.w[k] + dt * acc.w[k];
94+ const wNew = x.w[k] + dt * vNew;
95+ posErr = Math.max(posErr, Math.abs(out.pos[2 * n + k] + out.pos[3 * n + k] - wNew));
96+ posScale = Math.max(posScale, Math.abs(dt * vNew));
97+ }
98+ results.push({
99+ name: c.name,
100+ accError: { u: err.u / scale.u, v: err.v / scale.v, w: err.w / scale.w },
101+ positionError: posErr / posScale,
102+ nodes: n,
103+ dt,
104+ });
105+ }
106+ return results;
107+}
108+
109+export function selfTestPassed(r: SelfTestCase): boolean {
110+ return r.accError.u < ACC_TOLERANCE && r.accError.v < ACC_TOLERANCE && r.accError.w < ACC_TOLERANCE && r.positionError < POSITION_TOLERANCE;
111+}
src/sim/cpuRef.tsadded+237−0View file
@@ -0,0 +1,237 @@
1+// Float64 reference implementation of the spatial discretization. It mirrors
2+// step.wgsl and is used for the time-step estimate, the tests, and the
3+// GPU-versus-CPU self-test. It is not used in the time loop.
4+
5+import { type Grid, type Material, insideCell, insideNode, nodeWeight } from "../params.ts";
6+
7+export interface Fields {
8+ u: Float64Array;
9+ v: Float64Array;
10+ w: Float64Array;
11+}
12+
13+export function fields(n: number): Fields {
14+ return { u: new Float64Array(n), v: new Float64Array(n), w: new Float64Array(n) };
15+}
16+
17+const G3 = 1 / Math.sqrt(3);
18+
19+/** Shape function derivatives of the bilinear cell at Gauss point gp (0..3), scaled by 2/h. */
20+function shapeDerivatives(gp: number, hx: number, hy: number, dNdx: Float64Array, dNdy: Float64Array): void {
21+ const xi = gp & 1 ? G3 : -G3;
22+ const eta = gp & 2 ? G3 : -G3;
23+ dNdx[0] = (-(1 - eta) / 4) * (2 / hx);
24+ dNdx[1] = ((1 - eta) / 4) * (2 / hx);
25+ dNdx[2] = ((1 + eta) / 4) * (2 / hx);
26+ dNdx[3] = (-(1 + eta) / 4) * (2 / hx);
27+ dNdy[0] = (-(1 - xi) / 4) * (2 / hy);
28+ dNdy[1] = (-(1 + xi) / 4) * (2 / hy);
29+ dNdy[2] = ((1 + xi) / 4) * (2 / hy);
30+ dNdy[3] = ((1 - xi) / 4) * (2 / hy);
31+}
32+
33+/**
34+ * Undamped accelerations (out.u, out.v, out.w) of every node for the
35+ * displacement fields f and the stress-free shape w0. Bending forces are the
36+ * gradient of the discrete Kirchhoff energy (natural free-edge conditions);
37+ * membrane forces come from bilinear cells with 2x2 Gauss integration and
38+ * the Marguerre form of the von Karman strains. Nodes outside the sheet get zero.
39+ */
40+export function accelerations(g: Grid, m: Material, f: Fields, out: Fields, w0?: Float64Array): void {
41+ const { nx, ny, hx, hy } = g;
42+ const { u, v, w } = f;
43+ const hx2 = hx * hx, hy2 = hy * hy;
44+ const nu = m.nu;
45+ const id = (i: number, j: number) => j * nx + i;
46+
47+ // A second difference exists where the node and both of its neighbors belong to the sheet.
48+ const hasA = (i: number, j: number) => insideNode(g, i - 1, j) && insideNode(g, i, j) && insideNode(g, i + 1, j);
49+ const hasB = (i: number, j: number) => insideNode(g, i, j - 1) && insideNode(g, i, j) && insideNode(g, i, j + 1);
50+ const dxx = (i: number, j: number) => (w[id(i + 1, j)] - 2 * w[id(i, j)] + w[id(i - 1, j)]) / hx2;
51+ const dyy = (i: number, j: number) => (w[id(i, j + 1)] - 2 * w[id(i, j)] + w[id(i, j - 1)]) / hy2;
52+ // Weighted bending moments at nodes; zero where the second difference is undefined.
53+ const Ma = (i: number, j: number): number => {
54+ if (!hasA(i, j)) return 0;
55+ return nodeWeight(g, i, j) * (dxx(i, j) + nu * (hasB(i, j) ? dyy(i, j) : 0));
56+ };
57+ const Mb = (i: number, j: number): number => {
58+ if (!hasB(i, j)) return 0;
59+ return nodeWeight(g, i, j) * (dyy(i, j) + nu * (hasA(i, j) ? dxx(i, j) : 0));
60+ };
61+ // Mixed second difference at a cell center; zero for cells outside the sheet.
62+ const Cc = (ci: number, cj: number): number => {
63+ if (!insideCell(g, ci, cj)) return 0;
64+ return (w[id(ci + 1, cj + 1)] - w[id(ci + 1, cj)] - w[id(ci, cj + 1)] + w[id(ci, cj)]) / (hx * hy);
65+ };
66+
67+ const q = (hx * hy) / 4; // area per Gauss point
68+ const U = new Float64Array(4), V = new Float64Array(4), Wn = new Float64Array(4), W0 = new Float64Array(4);
69+ const dNdx = new Float64Array(4), dNdy = new Float64Array(4);
70+
71+ // Adds the membrane force of cell (ci, cj) on its local node a to F = [Fu, Fv, Fw].
72+ const cellForce = (ci: number, cj: number, a: number, F: Float64Array): void => {
73+ if (!insideCell(g, ci, cj)) return;
74+ const nn = [id(ci, cj), id(ci + 1, cj), id(ci + 1, cj + 1), id(ci, cj + 1)];
75+ for (let b = 0; b < 4; b++) {
76+ U[b] = u[nn[b]]; V[b] = v[nn[b]]; Wn[b] = w[nn[b]]; W0[b] = w0 ? w0[nn[b]] : 0;
77+ }
78+ for (let gp = 0; gp < 4; gp++) {
79+ shapeDerivatives(gp, hx, hy, dNdx, dNdy);
80+ let ux = 0, uy = 0, vx = 0, vy = 0, wx = 0, wy = 0, w0x = 0, w0y = 0;
81+ for (let b = 0; b < 4; b++) {
82+ ux += dNdx[b] * U[b]; uy += dNdy[b] * U[b];
83+ vx += dNdx[b] * V[b]; vy += dNdy[b] * V[b];
84+ wx += dNdx[b] * (Wn[b] + W0[b]); wy += dNdy[b] * (Wn[b] + W0[b]);
85+ w0x += dNdx[b] * W0[b]; w0y += dNdy[b] * W0[b];
86+ }
87+ // Marguerre strains: total slope (w0 + w) minus the stress-free slope of w0.
88+ const exx = ux + 0.5 * (wx * wx - w0x * w0x);
89+ const eyy = vy + 0.5 * (wy * wy - w0y * w0y);
90+ const gxy = uy + vx + wx * wy - w0x * w0y;
91+ const Nxx = m.C * (exx + nu * eyy);
92+ const Nyy = m.C * (eyy + nu * exx);
93+ const Nxy = m.G * gxy;
94+ F[0] -= q * (Nxx * dNdx[a] + Nxy * dNdy[a]);
95+ F[1] -= q * (Nyy * dNdy[a] + Nxy * dNdx[a]);
96+ F[2] -= q * ((Nxx * wx + Nxy * wy) * dNdx[a] + (Nxy * wx + Nyy * wy) * dNdy[a]);
97+ }
98+ };
99+
100+ const F = new Float64Array(3);
101+ for (let j = 0; j < ny; j++) {
102+ for (let i = 0; i < nx; i++) {
103+ const k = id(i, j);
104+ const weight = nodeWeight(g, i, j);
105+ if (weight === 0) {
106+ out.u[k] = 0; out.v[k] = 0; out.w[k] = 0;
107+ continue;
108+ }
109+ const bend =
110+ (Ma(i - 1, j) - 2 * Ma(i, j) + Ma(i + 1, j)) / hx2 +
111+ (Mb(i, j - 1) - 2 * Mb(i, j) + Mb(i, j + 1)) / hy2 +
112+ (2 * (1 - nu) * (Cc(i, j) - Cc(i - 1, j) - Cc(i, j - 1) + Cc(i - 1, j - 1))) / (hx * hy);
113+ F[0] = 0; F[1] = 0; F[2] = -m.D * hx * hy * bend;
114+ cellForce(i, j, 0, F);
115+ cellForce(i - 1, j, 1, F);
116+ cellForce(i - 1, j - 1, 2, F);
117+ cellForce(i, j - 1, 3, F);
118+ const mass = m.rhoh * hx * hy * weight;
119+ out.u[k] = F[0] / mass;
120+ out.v[k] = F[1] / mass;
121+ out.w[k] = F[2] / mass;
122+ }
123+ }
124+}
125+
126+/**
127+ * Adds the loss terms to the accelerations `out` for the velocities `vel`:
128+ * -2 sigma0 v on every component, +2 sigma1 (graph Laplacian of w_t) on the
129+ * transverse component, and -2 sigmaIn v on the in-plane components.
130+ */
131+export function addDamping(g: Grid, sigma0: number, sigma1: number, sigmaIn: number, vel: Fields, out: Fields): void {
132+ const { nx, ny, hx, hy } = g;
133+ const id = (i: number, j: number) => j * nx + i;
134+ for (let j = 0; j < ny; j++) {
135+ for (let i = 0; i < nx; i++) {
136+ const k = id(i, j);
137+ if (nodeWeight(g, i, j) === 0) continue;
138+ const c = vel.w[k];
139+ let lap = 0;
140+ if (insideNode(g, i - 1, j)) lap += (vel.w[id(i - 1, j)] - c) / (hx * hx);
141+ if (insideNode(g, i + 1, j)) lap += (vel.w[id(i + 1, j)] - c) / (hx * hx);
142+ if (insideNode(g, i, j - 1)) lap += (vel.w[id(i, j - 1)] - c) / (hy * hy);
143+ if (insideNode(g, i, j + 1)) lap += (vel.w[id(i, j + 1)] - c) / (hy * hy);
144+ out.u[k] -= 2 * (sigma0 + sigmaIn) * vel.u[k];
145+ out.v[k] -= 2 * (sigma0 + sigmaIn) * vel.v[k];
146+ out.w[k] += -2 * sigma0 * c + 2 * sigma1 * lap;
147+ }
148+ }
149+}
150+
151+/** Discrete potential energy (bending + membrane) whose negative gradient is the force in `accelerations`. */
152+export function potentialEnergy(g: Grid, m: Material, f: Fields, w0?: Float64Array): number {
153+ const { nx, ny, hx, hy } = g;
154+ const { u, v, w } = f;
155+ const nu = m.nu;
156+ const id = (i: number, j: number) => j * nx + i;
157+ let Ub = 0;
158+ for (let j = 0; j < ny; j++) {
159+ for (let i = 0; i < nx; i++) {
160+ let a = 0, b = 0;
161+ if (insideNode(g, i - 1, j) && insideNode(g, i, j) && insideNode(g, i + 1, j)) a = (w[id(i + 1, j)] - 2 * w[id(i, j)] + w[id(i - 1, j)]) / (hx * hx);
162+ if (insideNode(g, i, j - 1) && insideNode(g, i, j) && insideNode(g, i, j + 1)) b = (w[id(i, j + 1)] - 2 * w[id(i, j)] + w[id(i, j - 1)]) / (hy * hy);
163+ Ub += nodeWeight(g, i, j) * (a * a + b * b + 2 * nu * a * b);
164+ }
165+ }
166+ for (let cj = 0; cj < ny - 1; cj++) {
167+ for (let ci = 0; ci < nx - 1; ci++) {
168+ if (!insideCell(g, ci, cj)) continue;
169+ const c = (w[id(ci + 1, cj + 1)] - w[id(ci + 1, cj)] - w[id(ci, cj + 1)] + w[id(ci, cj)]) / (hx * hy);
170+ Ub += 2 * (1 - nu) * c * c;
171+ }
172+ }
173+ Ub *= 0.5 * m.D * hx * hy;
174+
175+ let Um = 0;
176+ const q = (hx * hy) / 4;
177+ const dNdx = new Float64Array(4), dNdy = new Float64Array(4);
178+ for (let cj = 0; cj < ny - 1; cj++) {
179+ for (let ci = 0; ci < nx - 1; ci++) {
180+ if (!insideCell(g, ci, cj)) continue;
181+ const nn = [id(ci, cj), id(ci + 1, cj), id(ci + 1, cj + 1), id(ci, cj + 1)];
182+ for (let gp = 0; gp < 4; gp++) {
183+ shapeDerivatives(gp, hx, hy, dNdx, dNdy);
184+ let ux = 0, uy = 0, vx = 0, vy = 0, wx = 0, wy = 0, w0x = 0, w0y = 0;
185+ for (let b = 0; b < 4; b++) {
186+ const z0 = w0 ? w0[nn[b]] : 0;
187+ ux += dNdx[b] * u[nn[b]]; uy += dNdy[b] * u[nn[b]];
188+ vx += dNdx[b] * v[nn[b]]; vy += dNdy[b] * v[nn[b]];
189+ wx += dNdx[b] * (w[nn[b]] + z0); wy += dNdy[b] * (w[nn[b]] + z0);
190+ w0x += dNdx[b] * z0; w0y += dNdy[b] * z0;
191+ }
192+ const exx = ux + 0.5 * (wx * wx - w0x * w0x);
193+ const eyy = vy + 0.5 * (wy * wy - w0y * w0y);
194+ const gxy = uy + vx + wx * wy - w0x * w0y;
195+ Um += q * 0.5 * (m.C * (exx * exx + eyy * eyy + 2 * nu * exx * eyy) + m.G * gxy * gxy);
196+ }
197+ }
198+ }
199+ return Ub + Um;
200+}
201+
202+/**
203+ * Largest angular frequency of the linearized discrete system, by power
204+ * iteration on a small full grid with the same spacing. The extreme mode is a
205+ * grid-scale in-plane mode, so neither the sheet size nor its outline matters.
206+ */
207+export function maxAngularFrequency(g: Grid, m: Material, nSmall = 12, iters = 300): number {
208+ const gs: Grid = { nx: nSmall, ny: nSmall, hx: g.hx, hy: g.hy, n: nSmall * nSmall };
209+ const n = gs.n;
210+ const x = fields(n), y = fields(n);
211+ let seed = 12345;
212+ const rnd = () => {
213+ seed = (seed * 1664525 + 1013904223) >>> 0;
214+ return seed / 4294967296 - 0.5;
215+ };
216+ for (let k = 0; k < n; k++) {
217+ x.u[k] = rnd(); x.v[k] = rnd(); x.w[k] = rnd();
218+ }
219+ const eps = 1e-9; // keeps the nonlinear terms negligible
220+ let lambda = 0;
221+ for (let it = 0; it < iters; it++) {
222+ let norm = 0;
223+ for (let k = 0; k < n; k++) norm += x.u[k] ** 2 + x.v[k] ** 2 + x.w[k] ** 2;
224+ norm = Math.sqrt(norm);
225+ for (let k = 0; k < n; k++) {
226+ x.u[k] *= eps / norm; x.v[k] *= eps / norm; x.w[k] *= eps / norm;
227+ }
228+ accelerations(gs, m, x, y);
229+ let dot = 0;
230+ for (let k = 0; k < n; k++) dot += x.u[k] * y.u[k] + x.v[k] * y.v[k] + x.w[k] * y.w[k];
231+ lambda = -dot / (eps * eps); // Rayleigh quotient of -M^{-1} K
232+ for (let k = 0; k < n; k++) {
233+ x.u[k] = -y.u[k]; x.v[k] = -y.v[k]; x.w[k] = -y.w[k];
234+ }
235+ }
236+ return Math.sqrt(Math.max(lambda, 0));
237+}
src/sim/shaders/audio.wgsladded+43−0View file
@@ -0,0 +1,43 @@
1+// Retarded-time sum over the plate for one output sample (Rayleigh integral).
2+// lw[k] = (weight, delay in samples) for node k; delay is relative to the
3+// closest node so that the ring buffer only has to cover the spread.
4+// The ring buffer holds the transverse velocity at the sample instants; the
5+// acceleration is the first difference of the delayed velocity, which is the
6+// average acceleration over one sample interval (band-limited, and immune to
7+// the per-step rounding noise of the force evaluation).
8+
9+@group(0) @binding(0) var<uniform> P: Params;
10+@group(0) @binding(1) var<uniform> SP: StepParams;
11+@group(0) @binding(2) var<storage, read> hist: array<f32>;
12+@group(0) @binding(3) var<storage, read> lw: array<vec2f>;
13+@group(0) @binding(4) var<storage, read_write> audioOut: array<f32>;
14+
15+var<workgroup> sh: array<f32, 256>;
16+
17+@compute @workgroup_size(256)
18+fn main(@builtin(local_invocation_id) lid: vec3u) {
19+ let t = lid.x;
20+ let H = i32(P.histLen);
21+ let slot = SP.audioSlot;
22+ var acc = 0.0;
23+ for (var k = t; k < P.n; k += 256u) {
24+ let d = lw[k].y;
25+ let k0 = i32(floor(d));
26+ let f = d - f32(k0);
27+ let s0 = ((slot - k0) % H + H) % H;
28+ let s1 = ((slot - k0 - 1) % H + H) % H;
29+ let s2 = ((slot - k0 - 2) % H + H) % H;
30+ let v0 = hist[u32(s0) * P.n + k];
31+ let v1 = hist[u32(s1) * P.n + k];
32+ let v2 = hist[u32(s2) * P.n + k];
33+ let a = (mix(v0, v1, f) - mix(v1, v2, f)) * P.fs;
34+ acc += lw[k].x * a;
35+ }
36+ sh[t] = acc;
37+ workgroupBarrier();
38+ for (var stride = 128u; stride > 0u; stride >>= 1u) {
39+ if (t < stride) { sh[t] += sh[t + stride]; }
40+ workgroupBarrier();
41+ }
42+ if (t == 0u) { audioOut[u32(SP.audioOut)] = sh[0]; }
43+}
src/sim/shaders/common.wgsladded+64−0View file
@@ -0,0 +1,64 @@
1+// Shared declarations, prepended to every kernel.
2+
3+struct Params {
4+ nx: u32,
5+ ny: u32,
6+ n: u32,
7+ histLen: u32,
8+ hx: f32,
9+ hy: f32,
10+ dt: f32,
11+ nu: f32,
12+ D: f32, // bending stiffness [N m]
13+ C: f32, // E h / (1 - nu^2) [N/m]
14+ G: f32, // E h / (2 (1 + nu)) [N/m]
15+ rhoh: f32, // mass per unit area [kg/m^2]
16+ pushX: f32, // center of the pushed patch [m]
17+ pushY: f32,
18+ pushR2: f32, // squared patch radius [m^2]
19+ mass: f32, // total mass [kg]
20+ Ix: f32, // sum m y^2
21+ Iy: f32, // sum m x^2
22+ Iz: f32, // sum m (x^2 + y^2)
23+ fs: f32, // audio sample rate [Hz]
24+ strikeX: f32, // drumstick contact point [m]
25+ strikeY: f32,
26+ strikeR2: f32, // squared contact radius [m^2]
27+ patchMass: f32, // mass of the contact patch [kg]
28+ strikeNode: u32, // sheet node nearest the contact point
29+ K: f32, // Hertz contact stiffness [N/m^alpha]
30+ alpha: f32, // Hertz exponent
31+ stickMass: f32, // [kg]
32+ patchI0: u32, // index box enclosing the contact patch (inclusive)
33+ patchI1: u32,
34+ patchJ0: u32,
35+ patchJ1: u32,
36+};
37+
38+// Per-step parameters, addressed with a dynamic uniform offset.
39+struct StepParams {
40+ step: u32,
41+ flags: u32, // bit 0: push constraint active; bit 1: record its reaction; bit 2: apply the recorded reaction; bit 3: drumstick active
42+ audioSlot: i32, // ring-buffer slot to record the velocity into (-1: none)
43+ audioOut: i32, // output sample index within the batch
44+ pushW: f32, // prescribed patch displacement [m]
45+ pushV: f32, // prescribed patch velocity [m/s]
46+ sigma0: f32,
47+ sigma1: f32,
48+ releaseGain: f32, // fraction of the recorded reaction still applied while the finger lets go
49+ sigmaIn: f32, // additional frequency-independent loss of the in-plane motion
50+ parity: u32, // which of the two stick-state slots holds the current state
51+ pad: u32,
52+};
53+
54+fn nodeX(i: i32) -> f32 { return (f32(i) - 0.5 * f32(P.nx - 1u)) * P.hx; }
55+fn nodeY(j: i32) -> f32 { return (f32(j) - 0.5 * f32(P.ny - 1u)) * P.hy; }
56+
57+// Quadrature and lumped-mass weight of every node: the fraction of its four
58+// cells that belong to the sheet (trapezoidal weights on a full rectangle),
59+// zero for nodes outside the sheet. Written once by the host.
60+@group(0) @binding(9) var<storage, read> nw: array<f32>;
61+fn nodeWeight(i: i32, j: i32) -> f32 {
62+ if (i < 0 || i >= i32(P.nx) || j < 0 || j >= i32(P.ny)) { return 0.0; }
63+ return nw[u32(j) * P.nx + u32(i)];
64+}
src/sim/shaders/extract.wgsladded+12−0View file
@@ -0,0 +1,12 @@
1+// Copies the transverse displacement into a plain f32 array (display / movie frame).
2+
3+@group(0) @binding(0) var<uniform> P: Params;
4+@group(0) @binding(1) var<storage, read> pos: array<f32>;
5+@group(0) @binding(2) var<storage, read_write> out: array<f32>;
6+
7+@compute @workgroup_size(256)
8+fn main(@builtin(global_invocation_id) gid: vec3u) {
9+ let k = gid.x;
10+ if (k >= P.n) { return; }
11+ out[k] = pos[2u * P.n + k];
12+}
src/sim/shaders/rigid.wgsladded+97−0View file
@@ -0,0 +1,97 @@
1+// Rigid-body projection: (1) mass-weighted sums of the six rigid modes for
2+// displacement and velocity, kinetic energy and max |w| in one workgroup;
3+// (2) subtraction of the rigid components from every node.
4+
5+@group(0) @binding(0) var<uniform> P: Params;
6+@group(0) @binding(1) var<storage, read_write> pos: array<f32>; // [u | v | w | wLo]
7+@group(0) @binding(2) var<storage, read_write> vel: array<f32>; // [du | dv | dw]
8+@group(0) @binding(3) var<storage, read_write> stats: array<f32>; // 0..11 coefficients, 12 KE, 13 max|w|
9+
10+const NS = 14u;
11+var<workgroup> sh: array<array<f32, NS>, 256>;
12+
13+@compute @workgroup_size(256)
14+fn sums(@builtin(local_invocation_id) lid: vec3u) {
15+ let t = lid.x;
16+ var s: array<f32, NS>;
17+ for (var q = 0u; q < NS; q++) { s[q] = 0.0; }
18+ for (var k = t; k < P.n; k += 256u) {
19+ let i = i32(k % P.nx); let j = i32(k / P.nx);
20+ let x = nodeX(i); let y = nodeY(j);
21+ let m = P.rhoh * P.hx * P.hy * nodeWeight(i, j);
22+ let n = P.n;
23+ let p = vec3f(pos[k], pos[n + k], pos[2u * n + k]);
24+ let v = vec3f(vel[k], vel[n + k], vel[2u * n + k]);
25+ s[0] += m * p.x; s[1] += m * p.y; s[2] += m * p.z;
26+ s[3] += m * (x * p.y - y * p.x); s[4] += m * y * p.z; s[5] += m * (-x * p.z);
27+ s[6] += m * v.x; s[7] += m * v.y; s[8] += m * v.z;
28+ s[9] += m * (x * v.y - y * v.x); s[10] += m * y * v.z; s[11] += m * (-x * v.z);
29+ s[12] += 0.5 * m * dot(v, v);
30+ s[13] = max(s[13], abs(p.z));
31+ }
32+ for (var q = 0u; q < NS; q++) { sh[t][q] = s[q]; }
33+ workgroupBarrier();
34+ for (var stride = 128u; stride > 0u; stride >>= 1u) {
35+ if (t < stride) {
36+ for (var q = 0u; q < 13u; q++) { sh[t][q] += sh[t + stride][q]; }
37+ sh[t][13] = max(sh[t][13], sh[t + stride][13]);
38+ }
39+ workgroupBarrier();
40+ }
41+ if (t == 0u) {
42+ let inv = array<f32, 6>(1.0 / P.mass, 1.0 / P.mass, 1.0 / P.mass, 1.0 / P.Iz, 1.0 / P.Ix, 1.0 / P.Iy);
43+ for (var q = 0u; q < 6u; q++) {
44+ stats[q] = sh[0][q] * inv[q];
45+ stats[q + 6u] = sh[0][q + 6u] * inv[q];
46+ }
47+ // Kinetic energy of the non-rigid motion (the six rigid modes are mass-orthogonal).
48+ var keRigid = 0.0;
49+ for (var q = 6u; q < 12u; q++) { keRigid += 0.5 * sh[0][q] * sh[0][q] * inv[q - 6u]; }
50+ stats[12] = max(sh[0][12] - keRigid, 0.0);
51+ stats[13] = sh[0][13];
52+ }
53+}
54+
55+@compute @workgroup_size(256)
56+fn apply(@builtin(global_invocation_id) gid: vec3u) {
57+ let k = gid.x;
58+ if (k >= P.n) { return; }
59+ let i = i32(k % P.nx); let j = i32(k / P.nx);
60+ if (nodeWeight(i, j) == 0.0) { return; } // outside the sheet
61+ let x = nodeX(i); let y = nodeY(j);
62+ let n = P.n;
63+ pos[k] -= stats[0] - stats[3] * y;
64+ pos[n + k] -= stats[1] + stats[3] * x;
65+ pos[2u * n + k] -= stats[2] + stats[4] * y - stats[5] * x;
66+ vel[k] -= stats[6] - stats[9] * y;
67+ vel[n + k] -= stats[7] + stats[9] * x;
68+ vel[2u * n + k] -= stats[8] + stats[10] * y - stats[11] * x;
69+}
70+
71+// Removes the three rigid components (translation and two rotations) from a
72+// per-node acceleration field so that it is self-equilibrated. One workgroup.
73+@group(0) @binding(4) var<storage, read_write> field: array<f32>;
74+var<workgroup> sh3: array<vec3f, 256>;
75+
76+@compute @workgroup_size(256)
77+fn projectField(@builtin(local_invocation_id) lid: vec3u) {
78+ let t = lid.x;
79+ var s = vec3f(0.0);
80+ for (var k = t; k < P.n; k += 256u) {
81+ let i = i32(k % P.nx); let j = i32(k / P.nx);
82+ let m = P.rhoh * P.hx * P.hy * nodeWeight(i, j);
83+ let f = field[k];
84+ s += vec3f(m * f, m * nodeY(j) * f, -m * nodeX(i) * f);
85+ }
86+ sh3[t] = s;
87+ workgroupBarrier();
88+ for (var stride = 128u; stride > 0u; stride >>= 1u) {
89+ if (t < stride) { sh3[t] += sh3[t + stride]; }
90+ workgroupBarrier();
91+ }
92+ let c = vec3f(sh3[0].x / P.mass, sh3[0].y / P.Ix, sh3[0].z / P.Iy);
93+ for (var k = t; k < P.n; k += 256u) {
94+ let i = i32(k % P.nx); let j = i32(k / P.nx);
95+ if (nodeWeight(i, j) > 0.0) { field[k] -= c.x + c.y * nodeY(j) - c.z * nodeX(i); }
96+ }
97+}
src/sim/shaders/step.wgsladded+267−0View file
@@ -0,0 +1,267 @@
1+// One explicit time step (symplectic Euler) of the von Karman plate.
2+//
3+// State layout (structure of arrays, n = nx * ny nodes):
4+// pos: [u | v | w | wLo] wLo carries the rounding error of the w update
5+// vel: [du | dv | dw]
6+// Each 16x16 workgroup loads its tile plus a 2-node halo of w and of the node
7+// weights into shared memory once, computes the membrane forces of the cells
8+// around the tile, and then evaluates the stencils from shared memory. Nodes
9+// outside the sheet (weight 0) are skipped: a second difference exists only
10+// where the node and both neighbors belong to the sheet, and a cell
11+// contributes only when its four corners do. That reproduces the free edge of
12+// the rectangle and gives a staircase free edge for an ellipse.
13+
14+@group(0) @binding(0) var<uniform> P: Params;
15+@group(0) @binding(1) var<uniform> SP: StepParams;
16+@group(0) @binding(2) var<storage, read> posIn: array<f32>;
17+@group(0) @binding(3) var<storage, read> velIn: array<f32>;
18+@group(0) @binding(4) var<storage, read_write> posOut: array<f32>;
19+@group(0) @binding(5) var<storage, read_write> velOut: array<f32>;
20+@group(0) @binding(6) var<storage, read_write> hist: array<f32>; // [slot * n + node]: transverse velocity at the audio sample instants
21+// Constraint reaction (acceleration of w) per node, followed by two 4-float
22+// slots of drumstick state (position, velocity, force, max force).
23+@group(0) @binding(7) var<storage, read_write> react: array<f32>;
24+@group(0) @binding(8) var<storage, read> w0: array<f32>; // stress-free initial shape
25+
26+const TILE = 16;
27+const HALO = 2;
28+const TW = 20; // TILE + 2 * HALO
29+var<workgroup> tw: array<f32, 400>;
30+// Node weights of the same tile in quarters (0..4), packed four to a word.
31+var<workgroup> tn: array<u32, 100>;
32+// Membrane force contributions of the 17x17 cells around the tile, 4 nodes x 3 components each.
33+const CT = 17;
34+var<workgroup> cellF: array<f32, 3468>;
35+
36+// Tile origin in global node coordinates (set once per invocation).
37+var<private> ox: i32;
38+var<private> oy: i32;
39+var<private> ihx2: f32;
40+var<private> ihy2: f32;
41+var<private> ihxy: f32;
42+
43+fn id(i: i32, j: i32) -> u32 { return u32(j) * P.nx + u32(i); }
44+fn lidx(i: i32, j: i32) -> u32 { return u32(j - oy) * u32(TW) + u32(i - ox); }
45+// Tile lookups with global indices; the tile covers [ox, ox + TW) x [oy, oy + TW).
46+fn W(i: i32, j: i32) -> f32 { return tw[lidx(i, j)]; }
47+fn NWq(i: i32, j: i32) -> u32 { let l = lidx(i, j); return (tn[l >> 2u] >> ((l & 3u) * 8u)) & 0xFFu; }
48+fn NW(i: i32, j: i32) -> f32 { return 0.25 * f32(NWq(i, j)); }
49+fn INS(i: i32, j: i32) -> bool { return NWq(i, j) != 0u; }
50+fn insCell(ci: i32, cj: i32) -> bool { return INS(ci, cj) && INS(ci + 1, cj) && INS(ci + 1, cj + 1) && INS(ci, cj + 1); }
51+fn U(i: i32, j: i32) -> f32 { return posIn[id(i, j)]; }
52+fn V(i: i32, j: i32) -> f32 { return posIn[P.n + id(i, j)]; }
53+fn DW(i: i32, j: i32) -> f32 { return velIn[2u * P.n + id(i, j)]; }
54+fn W0(i: i32, j: i32) -> f32 { return w0[id(i, j)]; }
55+
56+// Graph Laplacian of the transverse velocity at node (i, j), over neighbors in the sheet.
57+fn lapDW(i: i32, j: i32) -> f32 {
58+ let c = DW(i, j);
59+ var lap = 0.0;
60+ if (INS(i - 1, j)) { lap += (DW(i - 1, j) - c) * ihx2; }
61+ if (INS(i + 1, j)) { lap += (DW(i + 1, j) - c) * ihx2; }
62+ if (INS(i, j - 1)) { lap += (DW(i, j - 1) - c) * ihy2; }
63+ if (INS(i, j + 1)) { lap += (DW(i, j + 1) - c) * ihy2; }
64+ return lap;
65+}
66+
67+// Weighted bending moments at nodes; zero where the second difference is undefined.
68+fn Ma(i: i32, j: i32) -> f32 {
69+ if (!(INS(i - 1, j) && INS(i, j) && INS(i + 1, j))) { return 0.0; }
70+ let a = (W(i + 1, j) - 2.0 * W(i, j) + W(i - 1, j)) * ihx2;
71+ var b = 0.0;
72+ if (INS(i, j - 1) && INS(i, j + 1)) { b = (W(i, j + 1) - 2.0 * W(i, j) + W(i, j - 1)) * ihy2; }
73+ return NW(i, j) * (a + P.nu * b);
74+}
75+fn Mb(i: i32, j: i32) -> f32 {
76+ if (!(INS(i, j - 1) && INS(i, j) && INS(i, j + 1))) { return 0.0; }
77+ let b = (W(i, j + 1) - 2.0 * W(i, j) + W(i, j - 1)) * ihy2;
78+ var a = 0.0;
79+ if (INS(i - 1, j) && INS(i + 1, j)) { a = (W(i + 1, j) - 2.0 * W(i, j) + W(i - 1, j)) * ihx2; }
80+ return NW(i, j) * (b + P.nu * a);
81+}
82+// Mixed second difference at a cell center; zero outside the sheet.
83+fn Cc(ci: i32, cj: i32) -> f32 {
84+ if (!insCell(ci, cj)) { return 0.0; }
85+ return (W(ci + 1, cj + 1) - W(ci + 1, cj) - W(ci, cj + 1) + W(ci, cj)) * ihxy;
86+}
87+
88+// Membrane forces of the bilinear cell (ci, cj) on its four nodes (2x2 Gauss),
89+// written to cellF at local cell index lc.
90+fn cellForces(ci: i32, cj: i32, lc: u32) {
91+ let base = lc * 12u;
92+ if (!insCell(ci, cj)) {
93+ for (var q = 0u; q < 12u; q++) { cellF[base + q] = 0.0; }
94+ return;
95+ }
96+ let Un = vec4f(U(ci, cj), U(ci + 1, cj), U(ci + 1, cj + 1), U(ci, cj + 1));
97+ let Vn = vec4f(V(ci, cj), V(ci + 1, cj), V(ci + 1, cj + 1), V(ci, cj + 1));
98+ let W0n = vec4f(W0(ci, cj), W0(ci + 1, cj), W0(ci + 1, cj + 1), W0(ci, cj + 1));
99+ let Wn = vec4f(W(ci, cj), W(ci + 1, cj), W(ci + 1, cj + 1), W(ci, cj + 1)) + W0n; // total shape
100+ let q = 0.25 * P.hx * P.hy;
101+ let ihx = 1.0 / P.hx; let ihy = 1.0 / P.hy;
102+ let g3 = 0.57735026918962576;
103+ var Fu = vec4f(0.0); var Fv = vec4f(0.0); var Fw = vec4f(0.0);
104+ for (var gp = 0; gp < 4; gp++) {
105+ var xi = -g3; if ((gp & 1) != 0) { xi = g3; }
106+ var eta = -g3; if ((gp & 2) != 0) { eta = g3; }
107+ let dNdx = vec4f(-(1.0 - eta), 1.0 - eta, 1.0 + eta, -(1.0 + eta)) * (0.5 * ihx);
108+ let dNdy = vec4f(-(1.0 - xi), -(1.0 + xi), 1.0 + xi, 1.0 - xi) * (0.5 * ihy);
109+ let ux = dot(dNdx, Un); let uy = dot(dNdy, Un);
110+ let vx = dot(dNdx, Vn); let vy = dot(dNdy, Vn);
111+ let wx = dot(dNdx, Wn); let wy = dot(dNdy, Wn);
112+ let w0x = dot(dNdx, W0n); let w0y = dot(dNdy, W0n);
113+ // Marguerre strains: total slope minus the stress-free slope.
114+ let exx = ux + 0.5 * (wx * wx - w0x * w0x);
115+ let eyy = vy + 0.5 * (wy * wy - w0y * w0y);
116+ let gxy = uy + vx + wx * wy - w0x * w0y;
117+ let Nxx = P.C * (exx + P.nu * eyy);
118+ let Nyy = P.C * (eyy + P.nu * exx);
119+ let Nxy = P.G * gxy;
120+ Fu -= q * (Nxx * dNdx + Nxy * dNdy);
121+ Fv -= q * (Nyy * dNdy + Nxy * dNdx);
122+ Fw -= q * ((Nxx * wx + Nxy * wy) * dNdx + (Nxy * wx + Nyy * wy) * dNdy);
123+ }
124+ for (var a = 0u; a < 4u; a++) {
125+ cellF[base + a * 3u] = Fu[a];
126+ cellF[base + a * 3u + 1u] = Fv[a];
127+ cellF[base + a * 3u + 2u] = Fw[a];
128+ }
129+}
130+
131+// True for the sheet nodes of the drumstick contact patch (same rule as contactPatch in params.ts).
132+fn inContactPatch(i: i32, j: i32) -> bool {
133+ if (nodeWeight(i, j) == 0.0) { return false; }
134+ let dxs = nodeX(i) - P.strikeX; let dys = nodeY(j) - P.strikeY;
135+ return dxs * dxs + dys * dys <= P.strikeR2 || id(i, j) == P.strikeNode;
136+}
137+
138+// Stored membrane force on local node `a` of cell (ci, cj).
139+fn cellForce(ci: i32, cj: i32, a: u32) -> vec3f {
140+ let lc = u32(cj - oy - 1) * u32(CT) + u32(ci - ox - 1);
141+ let base = lc * 12u + a * 3u;
142+ return vec3f(cellF[base], cellF[base + 1u], cellF[base + 2u]);
143+}
144+
145+@compute @workgroup_size(16, 16)
146+fn main(@builtin(global_invocation_id) gid: vec3u, @builtin(local_invocation_id) lid: vec3u, @builtin(workgroup_id) wid: vec3u) {
147+ let n = P.n;
148+ ox = i32(wid.x) * TILE - HALO;
149+ oy = i32(wid.y) * TILE - HALO;
150+ ihx2 = 1.0 / (P.hx * P.hx); ihy2 = 1.0 / (P.hy * P.hy); ihxy = 1.0 / (P.hx * P.hy);
151+ // Cooperative load of the w and weight tiles (400 entries each, 256 threads).
152+ let lt = lid.y * 16u + lid.x;
153+ for (var e = lt; e < 400u; e += 256u) {
154+ let gi = ox + i32(e % u32(TW)); let gj = oy + i32(e / u32(TW));
155+ var ww = 0.0;
156+ if (gi >= 0 && gi < i32(P.nx) && gj >= 0 && gj < i32(P.ny)) { ww = posIn[2u * n + id(gi, gj)]; }
157+ tw[e] = ww;
158+ }
159+ if (lt < 100u) {
160+ var word = 0u;
161+ for (var q = 0u; q < 4u; q++) {
162+ let e = lt * 4u + q;
163+ let gi = ox + i32(e % u32(TW)); let gj = oy + i32(e / u32(TW));
164+ var wq = 0u;
165+ if (gi >= 0 && gi < i32(P.nx) && gj >= 0 && gj < i32(P.ny)) { wq = u32(round(4.0 * nw[id(gi, gj)])); }
166+ word |= wq << (q * 8u);
167+ }
168+ tn[lt] = word;
169+ }
170+ workgroupBarrier();
171+ // Membrane forces of the 17x17 cells adjacent to the tile's nodes.
172+ for (var e = lt; e < 289u; e += 256u) {
173+ let lci = i32(e % u32(CT)); let lcj = i32(e / u32(CT));
174+ cellForces(ox + 1 + lci, oy + 1 + lcj, e);
175+ }
176+ workgroupBarrier();
177+
178+ let i = i32(gid.x); let j = i32(gid.y);
179+ if (i >= i32(P.nx) || j >= i32(P.ny)) { return; }
180+ let k = id(i, j);
181+ let weight = NW(i, j);
182+ if (weight == 0.0) {
183+ // Outside the sheet: the node stays at rest.
184+ posOut[k] = 0.0; posOut[n + k] = 0.0; posOut[2u * n + k] = 0.0; posOut[3u * n + k] = 0.0;
185+ velOut[k] = 0.0; velOut[n + k] = 0.0; velOut[2u * n + k] = 0.0;
186+ if ((SP.flags & 2u) != 0u) { react[k] = 0.0; }
187+ return;
188+ }
189+
190+ // Elastic forces: bending from the discrete Kirchhoff energy, membrane from the cells.
191+ let bend = (Ma(i - 1, j) - 2.0 * Ma(i, j) + Ma(i + 1, j)) * ihx2
192+ + (Mb(i, j - 1) - 2.0 * Mb(i, j) + Mb(i, j + 1)) * ihy2
193+ + 2.0 * (1.0 - P.nu) * (Cc(i, j) - Cc(i - 1, j) - Cc(i, j - 1) + Cc(i - 1, j - 1)) * ihxy;
194+ var F = vec3f(0.0, 0.0, -P.D * P.hx * P.hy * bend);
195+ F += cellForce(i, j, 0u);
196+ F += cellForce(i - 1, j, 1u);
197+ F += cellForce(i - 1, j - 1, 2u);
198+ F += cellForce(i, j - 1, 3u);
199+ let mass = P.rhoh * P.hx * P.hy * weight;
200+ var acc = F * (1.0 / mass);
201+
202+ // Losses: frequency-independent on all components, the graph Laplacian of the
203+ // transverse velocity (~k^2), and an extra frequency-independent in-plane loss.
204+ let vel = vec3f(velIn[k], velIn[n + k], velIn[2u * n + k]);
205+ acc -= 2.0 * SP.sigma0 * vel;
206+ acc.z += 2.0 * SP.sigma1 * lapDW(i, j);
207+ acc.x -= 2.0 * SP.sigmaIn * vel.x;
208+ acc.y -= 2.0 * SP.sigmaIn * vel.y;
209+
210+ // Finger letting go: the reaction recorded at the end of the hold, ramped to zero.
211+ if ((SP.flags & 4u) != 0u) { acc.z += react[k] * SP.releaseGain; }
212+
213+ // Drumstick: Hertz contact between a free mass and the contact patch (the
214+ // sheet nodes within the tip radius, plus the nearest one). The contact
215+ // coordinate is the mass-weighted mean displacement of the patch, and the
216+ // force is spread over the patch in proportion to the node masses, so the
217+ // pair is conjugate and the contact has a potential energy.
218+ if ((SP.flags & 8u) != 0u) {
219+ if (inContactPatch(i, j)) {
220+ let sIn = n + 4u * SP.parity;
221+ let sx = react[sIn]; let sv = react[sIn + 1u]; let sMax = react[sIn + 3u];
222+ var wc = 0.0;
223+ for (var pj = i32(P.patchJ0); pj <= i32(P.patchJ1); pj++) {
224+ for (var pi = i32(P.patchI0); pi <= i32(P.patchI1); pi++) {
225+ if (inContactPatch(pi, pj)) { wc += P.rhoh * P.hx * P.hy * nodeWeight(pi, pj) * posIn[2u * n + id(pi, pj)]; }
226+ }
227+ }
228+ wc /= P.patchMass;
229+ let pen = wc - sx; // the stick sits on the +z side and moves toward -z
230+ var Fc = 0.0;
231+ if (pen > 0.0) { Fc = P.K * pow(pen, P.alpha); }
232+ acc.z -= Fc / P.patchMass;
233+ if (k == P.strikeNode) {
234+ let sOut = n + 4u * (1u - SP.parity);
235+ let svn = sv + P.dt * Fc / P.stickMass;
236+ react[sOut] = sx + P.dt * svn;
237+ react[sOut + 1u] = svn;
238+ react[sOut + 2u] = Fc;
239+ react[sOut + 3u] = max(sMax, Fc);
240+ }
241+ }
242+ }
243+
244+ var vNew = vel + P.dt * acc;
245+ if (SP.audioSlot >= 0) { hist[u32(SP.audioSlot) * n + k] = vNew.z; }
246+ let pOld = vec3f(posIn[k], posIn[n + k], posIn[2u * n + k]);
247+ var pNew = pOld + P.dt * vNew;
248+ // Compensated update of w: the increment dt * v can be far below the
249+ // float32 resolution of w, so the rounding error is carried separately.
250+ let y = P.dt * vNew.z + posIn[3u * n + k];
251+ let t = pOld.z + y;
252+ var lo = (pOld.z - t) + y;
253+ pNew.z = t;
254+ if ((SP.flags & 1u) != 0u) {
255+ let dxp = nodeX(i) - P.pushX; let dyp = nodeY(j) - P.pushY;
256+ let inPatch = dxp * dxp + dyp * dyp <= P.pushR2;
257+ if ((SP.flags & 2u) != 0u) {
258+ // Acceleration the constraint has to supply to hold this node.
259+ var r = 0.0;
260+ if (inPatch) { r = (SP.pushV - vNew.z) / P.dt; }
261+ react[k] = r;
262+ }
263+ if (inPatch) { pNew.z = SP.pushW; vNew.z = SP.pushV; lo = 0.0; }
264+ }
265+ posOut[k] = pNew.x; posOut[n + k] = pNew.y; posOut[2u * n + k] = pNew.z; posOut[3u * n + k] = lo;
266+ velOut[k] = vNew.x; velOut[n + k] = vNew.y; velOut[2u * n + k] = vNew.z;
267+}
src/sim/simulation.tsadded+662−0View file
@@ -0,0 +1,662 @@
1+// WebGPU time-stepping driver for the von Karman plate.
2+//
3+// The simulation advances in batches of one movie frame (1/fps seconds). A
4+// batch encodes one compute dispatch per time step (with a per-step uniform
5+// addressed by dynamic offset), a rigid-body projection every RB_INTERVAL
6+// steps, one retarded-time listener sum per audio sample, and one movie-frame
7+// extraction at the end. Audio, frame, and statistics are read back after
8+// every batch.
9+
10+import {
11+ type PhysParams,
12+ type Grid,
13+ type Material,
14+ makeGrid,
15+ makeMaterial,
16+ initialShape,
17+ strikePoint,
18+ pushPoint,
19+ contactPatch,
20+ nodeX,
21+ nodeY,
22+ nodeWeight,
23+ nodeWeights,
24+ HIST_LEN,
25+ RB_INTERVAL,
26+ RB_INTERVAL_HELD,
27+ RHO_AIR,
28+ SPEED_OF_SOUND,
29+} from "../params.ts";
30+import { maxAngularFrequency } from "./cpuRef.ts";
31+import commonWgsl from "./shaders/common.wgsl?raw";
32+import stepWgsl from "./shaders/step.wgsl?raw";
33+import rigidWgsl from "./shaders/rigid.wgsl?raw";
34+import audioWgsl from "./shaders/audio.wgsl?raw";
35+import extractWgsl from "./shaders/extract.wgsl?raw";
36+
37+export type Phase = "push" | "hold" | "recording" | "done";
38+
39+export interface StickState {
40+ x: number;
41+ v: number;
42+ force: number;
43+ maxForce: number;
44+}
45+
46+export interface SimStatus {
47+ phase: Phase;
48+ /** Simulation time relative to the release or the impact [s]. */
49+ t: number;
50+ totalSteps: number;
51+ ke: number;
52+ wmax: number;
53+ samples: number;
54+ frames: number;
55+ peakPressure: number;
56+ delayClamped: boolean;
57+ /** Drumstick state (strike excitation only). */
58+ stick: StickState | null;
59+}
60+
61+const PARAMS_LEN = 32; // floats in the Params uniform
62+const STEP_LEN = 12; // floats in the StepParams uniform
63+const STATS_LEN = 16;
64+const STICK_LEN = 8;
65+
66+export class Simulation {
67+ readonly grid: Grid;
68+ readonly mat: Material;
69+ readonly dt: number;
70+ /** Time steps per audio sample. */
71+ readonly stepsPerSample: number;
72+ readonly samplesPerBatch: number;
73+ readonly stepsPerBatch: number;
74+ readonly omegaMax: number;
75+
76+ /** Listener pressure [Pa] at fs, from the release or the impact. */
77+ audio: Float32Array = new Float32Array(0);
78+ samples = 0;
79+ /** Movie frames of w [m]; frame k is at t = k / fps. */
80+ frames: Float32Array<ArrayBuffer>[] = [];
81+
82+ phase: Phase = "push";
83+ totalSteps = 0;
84+ ke = 0;
85+ wmax = 0;
86+ peakPressure = 0;
87+ delayClamped = false;
88+
89+ /** Steps since the start of the current clock (preparation or recording). */
90+ private stepCount = 0;
91+ private prepTime = 0;
92+ private released = false;
93+ private lastPrepFrame: Float32Array<ArrayBuffer> | null = null;
94+ private stick: StickState | null = null;
95+
96+ private device: GPUDevice;
97+ private params: PhysParams;
98+ private paramsBuf: GPUBuffer;
99+ private stepBuf: GPUBuffer;
100+ private stepStride: number;
101+ private stepData: ArrayBuffer;
102+ private pos: GPUBuffer[];
103+ private vel: GPUBuffer[];
104+ private hist: GPUBuffer;
105+ private react: GPUBuffer;
106+ /** Stress-free initial shape (f32 per node), also used by the renderer. */
107+ readonly w0: GPUBuffer;
108+ /** Node weights (f32 per node, 0 outside the sheet), also used by the renderer. */
109+ readonly nw: GPUBuffer;
110+ private lw: GPUBuffer;
111+ private audioOut: GPUBuffer;
112+ private frameBuf: GPUBuffer;
113+ private statsBuf: GPUBuffer;
114+ private staging: GPUBuffer[];
115+ private stagingIdx = 0;
116+ private stagingSize: number;
117+ /** Index of the buffer pair holding the current state. */
118+ private cur = 0;
119+ /** Plain f32 array of w for the live view. */
120+ readonly displayW: GPUBuffer;
121+
122+ private stepPipe: GPUComputePipeline;
123+ private sumsPipe: GPUComputePipeline;
124+ private applyPipe: GPUComputePipeline;
125+ private projectPipe: GPUComputePipeline;
126+ private audioPipe: GPUComputePipeline;
127+ private extractPipe: GPUComputePipeline;
128+ private stepBG: GPUBindGroup[];
129+ private rbBG: GPUBindGroup[];
130+ private projectBG: GPUBindGroup;
131+ private audioBG: GPUBindGroup;
132+ private extractBG: GPUBindGroup[];
133+ private displayBG: GPUBindGroup[];
134+
135+ constructor(device: GPUDevice, params: PhysParams) {
136+ this.device = device;
137+ this.params = { ...params };
138+ const g = makeGrid(params);
139+ const m = makeMaterial(params);
140+ this.grid = g;
141+ this.mat = m;
142+
143+ // Stable time step, rounded so that an integer number of steps make one audio sample.
144+ this.omegaMax = maxAngularFrequency(g, m);
145+ const dtMax = (0.85 * 2) / this.omegaMax;
146+ this.stepsPerSample = Math.ceil(1 / (params.fs * dtMax));
147+ this.dt = 1 / (params.fs * this.stepsPerSample);
148+ this.samplesPerBatch = Math.max(1, Math.round(params.fs / params.fps));
149+ this.stepsPerBatch = this.samplesPerBatch * this.stepsPerSample;
150+
151+ const n = g.n;
152+ const f32 = (size: number, usage = GPUBufferUsage.STORAGE | GPUBufferUsage.COPY_SRC | GPUBufferUsage.COPY_DST) =>
153+ device.createBuffer({ size: Math.max(16, size * 4), usage });
154+
155+ this.pos = [f32(4 * n), f32(4 * n)]; // [u | v | w | wLo]
156+ this.vel = [f32(3 * n), f32(3 * n)]; // [du | dv | dw]
157+ this.hist = f32(HIST_LEN * n);
158+ this.react = f32(n + STICK_LEN);
159+ this.w0 = f32(n);
160+ device.queue.writeBuffer(this.w0, 0, initialShape(params, g));
161+ this.nw = f32(n);
162+ device.queue.writeBuffer(this.nw, 0, nodeWeights(g));
163+ this.lw = f32(2 * n);
164+ this.audioOut = f32(this.samplesPerBatch);
165+ this.frameBuf = f32(n);
166+ this.displayW = f32(n);
167+ this.statsBuf = f32(STATS_LEN);
168+ this.stagingSize = (this.samplesPerBatch + n + STATS_LEN + STICK_LEN) * 4;
169+ this.staging = [0, 1].map(() =>
170+ device.createBuffer({ size: this.stagingSize, usage: GPUBufferUsage.MAP_READ | GPUBufferUsage.COPY_DST }),
171+ );
172+
173+ // Mass and moments of inertia of the sheet.
174+ let mass = 0, Ix = 0, Iy = 0;
175+ for (let j = 0; j < g.ny; j++) {
176+ for (let i = 0; i < g.nx; i++) {
177+ const mk = m.rhoh * g.hx * g.hy * nodeWeight(g, i, j);
178+ mass += mk;
179+ Ix += mk * nodeY(g, j) ** 2;
180+ Iy += mk * nodeX(g, i) ** 2;
181+ }
182+ }
183+ // Drumstick contact patch and initial stick state (position 0, moving toward the sheet).
184+ const cp = contactPatch(params, g);
185+ const stick = new Float32Array(STICK_LEN);
186+ stick[1] = -params.strikeVelocity; stick[5] = -params.strikeVelocity;
187+ device.queue.writeBuffer(this.react, n * 4, stick);
188+
189+ const [px, py] = pushPoint(params);
190+ const [sx, sy] = strikePoint(params);
191+ const pf = new Float32Array(PARAMS_LEN);
192+ const pu = new Uint32Array(pf.buffer);
193+ pu[0] = g.nx; pu[1] = g.ny; pu[2] = n; pu[3] = HIST_LEN;
194+ pf[4] = g.hx; pf[5] = g.hy; pf[6] = this.dt; pf[7] = m.nu;
195+ pf[8] = m.D; pf[9] = m.C; pf[10] = m.G; pf[11] = m.rhoh;
196+ pf[12] = px; pf[13] = py; pf[14] = params.fingerRadius ** 2; pf[15] = mass;
197+ pf[16] = Ix; pf[17] = Iy; pf[18] = Ix + Iy; pf[19] = params.fs;
198+ pf[20] = sx; pf[21] = sy; pf[22] = params.contactRadius ** 2; pf[23] = cp.mass;
199+ pu[24] = cp.nearest; pf[25] = params.contactStiffness; pf[26] = params.contactExponent; pf[27] = params.stickMass;
200+ pu[28] = cp.box[0]; pu[29] = cp.box[1]; pu[30] = cp.box[2]; pu[31] = cp.box[3];
201+ this.paramsBuf = device.createBuffer({ size: pf.byteLength, usage: GPUBufferUsage.UNIFORM | GPUBufferUsage.COPY_DST });
202+ device.queue.writeBuffer(this.paramsBuf, 0, pf);
203+
204+ // Per-step parameters, one struct per step at a dynamic offset.
205+ this.stepStride = Math.max(STEP_LEN * 4, device.limits.minUniformBufferOffsetAlignment);
206+ this.stepData = new ArrayBuffer(this.stepStride * this.stepsPerBatch);
207+ this.stepBuf = device.createBuffer({ size: this.stepData.byteLength, usage: GPUBufferUsage.UNIFORM | GPUBufferUsage.COPY_DST });
208+
209+ // Pipelines.
210+ const shader = (code: string) => device.createShaderModule({ code: commonWgsl + code });
211+ const entry = (binding: number, type: GPUBufferBindingType, hasDynamicOffset = false): GPUBindGroupLayoutEntry => ({
212+ binding,
213+ visibility: GPUShaderStage.COMPUTE,
214+ buffer: { type, hasDynamicOffset },
215+ });
216+ const pipeline = (module: GPUShaderModule, entryPoint: string, layout: GPUBindGroupLayout) =>
217+ device.createComputePipeline({ layout: device.createPipelineLayout({ bindGroupLayouts: [layout] }), compute: { module, entryPoint } });
218+
219+ const stepLayout = device.createBindGroupLayout({
220+ entries: [
221+ entry(0, "uniform"),
222+ entry(1, "uniform", true),
223+ entry(2, "read-only-storage"),
224+ entry(3, "read-only-storage"),
225+ entry(4, "storage"),
226+ entry(5, "storage"),
227+ entry(6, "storage"),
228+ entry(7, "storage"),
229+ entry(8, "read-only-storage"),
230+ entry(9, "read-only-storage"),
231+ ],
232+ });
233+ this.stepPipe = pipeline(shader(stepWgsl), "main", stepLayout);
234+
235+ const rigidModule = shader(rigidWgsl);
236+ const rbLayout = device.createBindGroupLayout({
237+ entries: [entry(0, "uniform"), entry(1, "storage"), entry(2, "storage"), entry(3, "storage"), entry(9, "read-only-storage")],
238+ });
239+ this.sumsPipe = pipeline(rigidModule, "sums", rbLayout);
240+ this.applyPipe = pipeline(rigidModule, "apply", rbLayout);
241+ const projectLayout = device.createBindGroupLayout({
242+ entries: [entry(0, "uniform"), entry(4, "storage"), entry(9, "read-only-storage")],
243+ });
244+ this.projectPipe = pipeline(rigidModule, "projectField", projectLayout);
245+
246+ const audioLayout = device.createBindGroupLayout({
247+ entries: [entry(0, "uniform"), entry(1, "uniform", true), entry(2, "read-only-storage"), entry(3, "read-only-storage"), entry(4, "storage")],
248+ });
249+ this.audioPipe = pipeline(shader(audioWgsl), "main", audioLayout);
250+
251+ const extractLayout = device.createBindGroupLayout({
252+ entries: [entry(0, "uniform"), entry(1, "read-only-storage"), entry(2, "storage")],
253+ });
254+ this.extractPipe = pipeline(shader(extractWgsl), "main", extractLayout);
255+
256+ // Bind groups. The state is double-buffered: bind group p reads pair p and writes pair 1 - p.
257+ const sp = { buffer: this.stepBuf, size: STEP_LEN * 4 };
258+ this.stepBG = [0, 1].map((p) =>
259+ device.createBindGroup({
260+ layout: stepLayout,
261+ entries: [
262+ { binding: 0, resource: { buffer: this.paramsBuf } },
263+ { binding: 1, resource: sp },
264+ { binding: 2, resource: { buffer: this.pos[p] } },
265+ { binding: 3, resource: { buffer: this.vel[p] } },
266+ { binding: 4, resource: { buffer: this.pos[1 - p] } },
267+ { binding: 5, resource: { buffer: this.vel[1 - p] } },
268+ { binding: 6, resource: { buffer: this.hist } },
269+ { binding: 7, resource: { buffer: this.react } },
270+ { binding: 8, resource: { buffer: this.w0 } },
271+ { binding: 9, resource: { buffer: this.nw } },
272+ ],
273+ }),
274+ );
275+ this.rbBG = [0, 1].map((p) =>
276+ device.createBindGroup({
277+ layout: rbLayout,
278+ entries: [
279+ { binding: 0, resource: { buffer: this.paramsBuf } },
280+ { binding: 1, resource: { buffer: this.pos[p] } },
281+ { binding: 2, resource: { buffer: this.vel[p] } },
282+ { binding: 3, resource: { buffer: this.statsBuf } },
283+ { binding: 9, resource: { buffer: this.nw } },
284+ ],
285+ }),
286+ );
287+ this.projectBG = device.createBindGroup({
288+ layout: projectLayout,
289+ entries: [
290+ { binding: 0, resource: { buffer: this.paramsBuf } },
291+ { binding: 4, resource: { buffer: this.react } },
292+ { binding: 9, resource: { buffer: this.nw } },
293+ ],
294+ });
295+ this.audioBG = device.createBindGroup({
296+ layout: audioLayout,
297+ entries: [
298+ { binding: 0, resource: { buffer: this.paramsBuf } },
299+ { binding: 1, resource: sp },
300+ { binding: 2, resource: { buffer: this.hist } },
301+ { binding: 3, resource: { buffer: this.lw } },
302+ { binding: 4, resource: { buffer: this.audioOut } },
303+ ],
304+ });
305+ const extractBG = (p: number, out: GPUBuffer) =>
306+ device.createBindGroup({
307+ layout: extractLayout,
308+ entries: [
309+ { binding: 0, resource: { buffer: this.paramsBuf } },
310+ { binding: 1, resource: { buffer: this.pos[p] } },
311+ { binding: 2, resource: { buffer: out } },
312+ ],
313+ });
314+ this.extractBG = [0, 1].map((p) => extractBG(p, this.frameBuf));
315+ this.displayBG = [0, 1].map((p) => extractBG(p, this.displayW));
316+
317+ this.setListener(params.listenerDist, params.listenerAzimuth, params.listenerElevation);
318+ if (params.excitation === "strike") {
319+ // No preparation: the sheet is at rest in its natural shape when the stick arrives.
320+ this.phase = "recording";
321+ this.released = true;
322+ this.frames.push(new Float32Array(n));
323+ }
324+ }
325+
326+ /** Recomputes the per-node weights and delays of the listener sum. Takes effect on later samples. */
327+ setListener(dist: number, azimuthDeg: number, elevationDeg: number): void {
328+ const g = this.grid;
329+ const az = (azimuthDeg * Math.PI) / 180;
330+ const el = (elevationDeg * Math.PI) / 180;
331+ const L = [dist * Math.sin(az) * Math.cos(el), dist * Math.sin(el), dist * Math.cos(az) * Math.cos(el)];
332+ const data = new Float32Array(2 * g.n);
333+ let dmin = Infinity;
334+ for (let j = 0; j < g.ny; j++) {
335+ for (let i = 0; i < g.nx; i++) {
336+ const k = j * g.nx + i;
337+ const rx = L[0] - nodeX(g, i), ry = L[1] - nodeY(g, j), rz = L[2];
338+ const r = Math.sqrt(rx * rx + ry * ry + rz * rz);
339+ const dA = g.hx * g.hy * nodeWeight(g, i, j);
340+ data[2 * k] = (RHO_AIR * dA) / (2 * Math.PI * r);
341+ data[2 * k + 1] = (r / SPEED_OF_SOUND) * this.params.fs;
342+ dmin = Math.min(dmin, data[2 * k + 1]);
343+ }
344+ }
345+ // Only the spread of delays across the sheet is kept; the common delay is dropped.
346+ this.delayClamped = false;
347+ for (let k = 0; k < g.n; k++) {
348+ let d = data[2 * k + 1] - dmin;
349+ if (d > HIST_LEN - 3) {
350+ d = HIST_LEN - 3;
351+ this.delayClamped = true;
352+ }
353+ data[2 * k + 1] = d;
354+ }
355+ this.device.queue.writeBuffer(this.lw, 0, data);
356+ this.params.listenerDist = dist;
357+ this.params.listenerAzimuth = azimuthDeg;
358+ this.params.listenerElevation = elevationDeg;
359+ }
360+
361+ setDamping(sigma0: number, sigma1: number, sigmaIn: number): void {
362+ this.params.sigma0 = sigma0;
363+ this.params.sigma1 = sigma1;
364+ this.params.sigmaIn = sigmaIn;
365+ }
366+
367+ setDuration(duration: number): void {
368+ this.params.duration = duration;
369+ if (this.canContinue()) this.phase = "recording";
370+ }
371+
372+ /** True when a finished or stopped run can be extended. */
373+ canContinue(): boolean {
374+ return this.phase === "done" && this.released && this.samples < this.targetSamples();
375+ }
376+
377+ /** Ends the run; the state, audio, and movie are kept. */
378+ stop(): void {
379+ this.phase = "done";
380+ }
381+
382+ targetSamples(): number {
383+ return Math.round(this.params.duration * this.params.fs);
384+ }
385+
386+ status(): SimStatus {
387+ const t = this.phase === "push" || this.phase === "hold" ? this.prepTime - this.prepDuration() : this.stepCount * this.dt;
388+ return {
389+ phase: this.phase,
390+ t,
391+ totalSteps: this.totalSteps,
392+ ke: this.ke,
393+ wmax: this.wmax,
394+ samples: this.samples,
395+ frames: this.frames.length,
396+ peakPressure: this.peakPressure,
397+ delayClamped: this.delayClamped,
398+ stick: this.stick,
399+ };
400+ }
401+
402+ private prepDuration(): number {
403+ return this.params.rampTime + this.params.holdMax;
404+ }
405+
406+ /** Enqueues an extraction of the current w field into `displayW` (for the live view). */
407+ encodeDisplay(encoder: GPUCommandEncoder): void {
408+ const pass = encoder.beginComputePass();
409+ pass.setPipeline(this.extractPipe);
410+ pass.setBindGroup(0, this.displayBG[this.cur]);
411+ pass.dispatchWorkgroups(Math.ceil(this.grid.n / 256));
412+ pass.end();
413+ }
414+
415+ /** The recorded listener pressure so far [Pa]. */
416+ recordedAudio(): Float32Array {
417+ return this.audio.subarray(0, this.samples);
418+ }
419+
420+ private appendAudio(raw: Float32Array): void {
421+ const need = this.samples + raw.length;
422+ if (this.audio.length < need) {
423+ const a = new Float32Array(Math.max(need, 2 * this.audio.length, this.params.fs));
424+ a.set(this.audio.subarray(0, this.samples));
425+ this.audio = a;
426+ }
427+ this.audio.set(raw, this.samples);
428+ for (let k = 0; k < raw.length; k++) this.peakPressure = Math.max(this.peakPressure, Math.abs(raw[k]));
429+ this.samples = need;
430+ }
431+
432+ /** Writes the per-step uniform for step s of a batch. */
433+ private writeStepParams(s: number, fields: { flags: number; slot: number; outIdx: number; pushW: number; pushV: number; sigma0: number; releaseGain: number }): void {
434+ const o = (s * this.stepStride) / 4;
435+ const u32 = new Uint32Array(this.stepData);
436+ const i32 = new Int32Array(this.stepData);
437+ const f32 = new Float32Array(this.stepData);
438+ u32[o] = this.totalSteps + s;
439+ u32[o + 1] = fields.flags;
440+ i32[o + 2] = fields.slot;
441+ i32[o + 3] = fields.outIdx;
442+ f32[o + 4] = fields.pushW;
443+ f32[o + 5] = fields.pushV;
444+ f32[o + 6] = fields.sigma0;
445+ f32[o + 7] = this.params.sigma1;
446+ f32[o + 8] = fields.releaseGain;
447+ f32[o + 9] = this.params.sigmaIn;
448+ u32[o + 10] = (this.totalSteps + s) & 1;
449+ u32[o + 11] = 0;
450+ }
451+
452+ /** Runs one batch (one movie frame of time). Resolves when its results have been read back. */
453+ async runBatch(): Promise<void> {
454+ if (this.phase === "done") return;
455+ const p = this.params;
456+ const recording = this.phase === "recording";
457+ const stepsLeft = recording ? (this.targetSamples() - this.samples) * this.stepsPerSample : this.stepsPerBatch;
458+ const steps = Math.min(this.stepsPerBatch, Math.max(0, stepsLeft));
459+ if (steps === 0) {
460+ this.phase = "done";
461+ return;
462+ }
463+
464+ // Per-step parameters: the push ramp and settle damping before the release,
465+ // the release ramp or the drumstick and the audio bookkeeping afterwards.
466+ const audioSteps: number[] = [];
467+ for (let s = 0; s < steps; s++) {
468+ const local = this.stepCount + s;
469+ let flags = 0, pushW = 0, pushV = 0, sigma0 = p.sigma0, slot = -1, outIdx = -1, releaseGain = 0;
470+ if (!recording) {
471+ const t = this.prepTime + (s + 1) * this.dt; // state after this step
472+ const tau = Math.min(1, t / p.rampTime);
473+ const smooth = tau * tau * (3 - 2 * tau);
474+ const dsmooth = tau < 1 ? (6 * tau * (1 - tau)) / p.rampTime : 0;
475+ flags = 1 | (tau >= 1 ? 2 : 0); // record the reaction once the ramp is complete
476+ pushW = p.pushAmount * smooth;
477+ pushV = p.pushAmount * dsmooth;
478+ sigma0 = p.settleSigma0;
479+ } else {
480+ const t = local * this.dt;
481+ if (p.excitation === "strike") {
482+ flags = 8;
483+ } else if (p.releaseTime > 0 && t < p.releaseTime) {
484+ const tau = t / p.releaseTime;
485+ releaseGain = 1 - tau * tau * (3 - 2 * tau);
486+ flags = 4;
487+ }
488+ if (local % this.stepsPerSample === 0) {
489+ slot = (local / this.stepsPerSample) % HIST_LEN;
490+ outIdx = audioSteps.length;
491+ audioSteps.push(s);
492+ }
493+ }
494+ this.writeStepParams(s, { flags, slot, outIdx, pushW, pushV, sigma0, releaseGain });
495+ }
496+ this.device.queue.writeBuffer(this.stepBuf, 0, this.stepData, 0, steps * this.stepStride);
497+
498+ // Encode the batch in chunks so that the GPU starts while JS is still encoding.
499+ const CHUNK = 1024;
500+ let enc = this.device.createCommandEncoder();
501+ let pass = enc.beginComputePass();
502+ const wgX = Math.ceil(this.grid.nx / 16), wgY = Math.ceil(this.grid.ny / 16);
503+ const wgN = Math.ceil(this.grid.n / 256);
504+ let cur = this.cur;
505+ let audioPtr = 0;
506+ const rbInterval = recording ? RB_INTERVAL : RB_INTERVAL_HELD;
507+ let stepBound = false;
508+ for (let s = 0; s < steps; s++) {
509+ if (s > 0 && s % CHUNK === 0) {
510+ pass.end();
511+ this.device.queue.submit([enc.finish()]);
512+ enc = this.device.createCommandEncoder();
513+ pass = enc.beginComputePass();
514+ stepBound = false;
515+ }
516+ if (!stepBound) {
517+ pass.setPipeline(this.stepPipe);
518+ stepBound = true;
519+ }
520+ pass.setBindGroup(0, this.stepBG[cur], [s * this.stepStride]);
521+ pass.dispatchWorkgroups(wgX, wgY);
522+ cur = 1 - cur;
523+ const global = this.totalSteps + s + 1;
524+ if (global % rbInterval === 0 || s === steps - 1) {
525+ pass.setPipeline(this.sumsPipe);
526+ pass.setBindGroup(0, this.rbBG[cur]);
527+ pass.dispatchWorkgroups(1);
528+ pass.setPipeline(this.applyPipe);
529+ pass.setBindGroup(0, this.rbBG[cur]);
530+ pass.dispatchWorkgroups(wgN);
531+ stepBound = false;
532+ }
533+ if (audioPtr < audioSteps.length && audioSteps[audioPtr] === s) {
534+ pass.setPipeline(this.audioPipe);
535+ pass.setBindGroup(0, this.audioBG, [s * this.stepStride]);
536+ pass.dispatchWorkgroups(1);
537+ audioPtr++;
538+ stepBound = false;
539+ }
540+ }
541+ pass.setPipeline(this.extractPipe);
542+ pass.setBindGroup(0, this.extractBG[cur]);
543+ pass.dispatchWorkgroups(wgN);
544+ pass.end();
545+ const staging = this.staging[this.stagingIdx];
546+ this.stagingIdx = 1 - this.stagingIdx;
547+ const nAudio = audioSteps.length;
548+ const n = this.grid.n;
549+ enc.copyBufferToBuffer(this.audioOut, 0, staging, 0, Math.max(4, nAudio * 4));
550+ enc.copyBufferToBuffer(this.frameBuf, 0, staging, this.samplesPerBatch * 4, n * 4);
551+ enc.copyBufferToBuffer(this.statsBuf, 0, staging, (this.samplesPerBatch + n) * 4, STATS_LEN * 4);
552+ enc.copyBufferToBuffer(this.react, n * 4, staging, (this.samplesPerBatch + n + STATS_LEN) * 4, STICK_LEN * 4);
553+ this.device.queue.submit([enc.finish()]);
554+ this.cur = cur;
555+ this.totalSteps += steps;
556+ this.stepCount += steps;
557+ if (!recording) this.prepTime += steps * this.dt;
558+
559+ await staging.mapAsync(GPUMapMode.READ);
560+ const view = new Float32Array(staging.getMappedRange());
561+ const stats = view.subarray(this.samplesPerBatch + n, this.samplesPerBatch + n + STATS_LEN);
562+ this.ke = stats[12];
563+ this.wmax = stats[13];
564+ if (p.excitation === "strike") {
565+ // The slot written by the last step is the one the next step reads.
566+ const base = this.samplesPerBatch + n + STATS_LEN + 4 * (this.totalSteps & 1);
567+ this.stick = { x: view[base], v: view[base + 1], force: view[base + 2], maxForce: view[base + 3] };
568+ }
569+ const frame = view.slice(this.samplesPerBatch, this.samplesPerBatch + n);
570+ if (recording) {
571+ this.appendAudio(view.slice(0, nAudio));
572+ this.frames.push(frame);
573+ } else {
574+ this.lastPrepFrame = frame;
575+ }
576+ staging.unmap();
577+
578+ if (!recording) {
579+ if (this.phase === "push" && this.prepTime >= p.rampTime) this.phase = "hold";
580+ const settled = this.phase === "hold" && this.ke < p.settleKE;
581+ if (settled || this.prepTime >= this.prepDuration()) this.release();
582+ } else if (this.samples >= this.targetSamples()) {
583+ this.phase = "done";
584+ }
585+ }
586+
587+ private release(): void {
588+ this.phase = "recording";
589+ this.released = true;
590+ this.stepCount = 0;
591+ this.totalSteps = 0;
592+ // The finger force and the distributed holding reaction relax together:
593+ // remove the rigid-body components of the recorded reaction.
594+ const enc = this.device.createCommandEncoder();
595+ const pass = enc.beginComputePass();
596+ pass.setPipeline(this.projectPipe);
597+ pass.setBindGroup(0, this.projectBG);
598+ pass.dispatchWorkgroups(1);
599+ pass.end();
600+ this.device.queue.submit([enc.finish()]);
601+ // Frame 0 is the held, settled shape.
602+ if (this.lastPrepFrame) this.frames.push(this.lastPrepFrame);
603+ }
604+
605+ // ---------- Testing support ----------
606+
607+ /**
608+ * Advances the state by exactly one free step (no constraint, no drumstick,
609+ * no rigid-body projection, no audio) with the current loss parameters.
610+ */
611+ async stepOnce(): Promise<void> {
612+ this.writeStepParams(0, { flags: 0, slot: -1, outIdx: -1, pushW: 0, pushV: 0, sigma0: this.params.sigma0, releaseGain: 0 });
613+ this.device.queue.writeBuffer(this.stepBuf, 0, this.stepData, 0, this.stepStride);
614+ const enc = this.device.createCommandEncoder();
615+ const pass = enc.beginComputePass();
616+ pass.setPipeline(this.stepPipe);
617+ pass.setBindGroup(0, this.stepBG[this.cur], [0]);
618+ pass.dispatchWorkgroups(Math.ceil(this.grid.nx / 16), Math.ceil(this.grid.ny / 16));
619+ pass.end();
620+ this.device.queue.submit([enc.finish()]);
621+ this.cur = 1 - this.cur;
622+ this.totalSteps++;
623+ await this.device.queue.onSubmittedWorkDone();
624+ }
625+
626+ /** Overwrites the state: pos as [u | v | w | wLo], vel as [du | dv | dw]. */
627+ writeState(pos: Float32Array<ArrayBuffer>, vel: Float32Array<ArrayBuffer>): void {
628+ const n = this.grid.n;
629+ if (pos.length !== 4 * n || vel.length !== 3 * n) throw new Error("state arrays have the wrong length");
630+ this.device.queue.writeBuffer(this.pos[this.cur], 0, pos);
631+ this.device.queue.writeBuffer(this.vel[this.cur], 0, vel);
632+ }
633+
634+ /** Reads back the full state as [u | v | w | wLo] and [du | dv | dw]. */
635+ async readState(): Promise<{ pos: Float32Array; vel: Float32Array }> {
636+ const n = this.grid.n;
637+ const stg = this.device.createBuffer({ size: 7 * n * 4, usage: GPUBufferUsage.MAP_READ | GPUBufferUsage.COPY_DST });
638+ const enc = this.device.createCommandEncoder();
639+ enc.copyBufferToBuffer(this.pos[this.cur], 0, stg, 0, 4 * n * 4);
640+ enc.copyBufferToBuffer(this.vel[this.cur], 0, stg, 4 * n * 4, 3 * n * 4);
641+ this.device.queue.submit([enc.finish()]);
642+ await stg.mapAsync(GPUMapMode.READ);
643+ const all = new Float32Array(stg.getMappedRange()).slice();
644+ stg.unmap();
645+ stg.destroy();
646+ return { pos: all.subarray(0, 4 * n), vel: all.subarray(4 * n) };
647+ }
648+
649+ destroy(): void {
650+ for (const b of [...this.pos, ...this.vel, this.hist, this.react, this.w0, this.nw, this.lw, this.audioOut, this.frameBuf, this.displayW, this.statsBuf, this.paramsBuf, this.stepBuf, ...this.staging]) b.destroy();
651+ }
652+
653+ get fs(): number {
654+ return this.params.fs;
655+ }
656+ get fps(): number {
657+ return this.params.fps;
658+ }
659+ get physParams(): PhysParams {
660+ return this.params;
661+ }
662+}
src/style.cssadded+58−0View file
@@ -0,0 +1,58 @@
1+:root {
2+ color-scheme: dark;
3+ --bg: #15171a;
4+ --panel: #1c1e22;
5+ --border: #2c3036;
6+ --text: #dfe2e6;
7+ --muted: #9aa0a8;
8+ --accent: #f0b429;
9+ --accent2: #3b8ce6;
10+ font-family: system-ui, -apple-system, "Segoe UI", Roboto, sans-serif;
11+ font-size: 13px;
12+}
13+* { box-sizing: border-box; }
14+html, body { margin: 0; height: 100%; background: var(--bg); color: var(--text); }
15+#app { display: flex; height: 100vh; }
16+#sidebar {
17+ width: 330px; flex: none; overflow-y: auto; padding: 12px 14px; background: var(--panel);
18+ border-right: 1px solid var(--border);
19+}
20+#main { flex: 1; display: flex; flex-direction: column; min-width: 0; padding: 8px; gap: 6px; }
21+h1 { font-size: 17px; margin: 0 0 4px; }
22+h2 { font-size: 12px; text-transform: uppercase; letter-spacing: 0.06em; color: var(--muted); margin: 14px 0 4px; }
23+.sub { color: var(--muted); margin: 0 0 8px; line-height: 1.35; }
24+a { color: var(--accent2); }
25+details { border-top: 1px solid var(--border); padding: 4px 0; }
26+details summary { cursor: pointer; font-size: 12px; text-transform: uppercase; letter-spacing: 0.06em; color: var(--muted); padding: 4px 0; }
27+.field { display: grid; grid-template-columns: 1fr 96px 34px; align-items: center; gap: 6px; margin: 3px 0; }
28+.field label { color: var(--text); }
29+.field .unit { color: var(--muted); font-size: 11px; }
30+input[type=number], select {
31+ width: 100%; background: #0f1113; color: var(--text); border: 1px solid var(--border); border-radius: 4px; padding: 3px 5px; font: inherit;
32+}
33+input[type=number]:disabled { color: var(--muted); }
34+button {
35+ background: #2a2e34; color: var(--text); border: 1px solid #3a3f47; border-radius: 5px; padding: 5px 10px; font: inherit; cursor: pointer;
36+}
37+button:hover { background: #343941; }
38+button:disabled { opacity: 0.45; cursor: default; }
39+button.primary { background: var(--accent); color: #1a1500; border-color: var(--accent); font-weight: 600; }
40+button.primary:hover { background: #f8c453; }
41+#controls { display: flex; flex-wrap: wrap; gap: 6px; margin: 8px 0; }
42+#info { margin-top: 10px; border-top: 1px solid var(--border); padding-top: 8px; }
43+#info .row { display: flex; justify-content: space-between; padding: 1px 0; }
44+#info .row span:last-child { font-variant-numeric: tabular-nums; color: var(--text); }
45+#info .row span:first-child { color: var(--muted); }
46+.warn { color: var(--accent); }
47+.err { color: #ff6b6b; }
48+.toolbar { display: flex; flex-wrap: wrap; align-items: center; gap: 8px; padding: 4px 6px; background: var(--panel); border: 1px solid var(--border); border-radius: 6px; }
49+.toolbar label { display: inline-flex; align-items: center; gap: 4px; color: var(--muted); }
50+.toolbar input[type=number], .toolbar select { width: 82px; }
51+.toolbar input[type=range] { flex: 1; min-width: 120px; }
52+.toolbar .time { font-variant-numeric: tabular-nums; min-width: 72px; }
53+#viewWrap { position: relative; flex: 3; min-height: 200px; border: 1px solid var(--border); border-radius: 6px; overflow: hidden; }
54+#view { width: 100%; height: 100%; display: block; touch-action: none; }
55+#spec { flex: 2; min-height: 140px; width: 100%; border: 1px solid var(--border); border-radius: 6px; }
56+#nogpu[hidden] { display: none; }
57+#nogpu { position: absolute; inset: 0; display: flex; align-items: center; justify-content: center; padding: 24px; text-align: center; color: #ff8a8a; background: var(--bg); }
58+.seg button.on { background: var(--accent2); border-color: var(--accent2); color: white; }
src/ui/form.tsadded+322−0View file
@@ -0,0 +1,322 @@
1+// Parameter form: numeric fields in convenient units, a few selects, and the
2+// derived-quantity notes. Everything is stored in SI in `params`.
3+
4+import {
5+ type BowShape,
6+ type Excitation,
7+ type Grid,
8+ type PhysParams,
9+ type PlateShape,
10+ countInside,
11+ dispersionFrequency,
12+ makeGrid,
13+ makeMaterial,
14+ t60,
15+} from "../params.ts";
16+
17+export type SizeUnit = "in" | "mm";
18+
19+interface FieldSpec {
20+ key: keyof PhysParams;
21+ label: string;
22+ unit: string;
23+ scale: number; // form value = SI value / scale
24+ step?: number;
25+ live?: boolean; // may be changed while a simulation runs
26+}
27+
28+interface SelectSpec {
29+ key: keyof PhysParams;
30+ label: string;
31+ options: [string, string][];
32+ before?: keyof PhysParams; // insert before this field (default: at the start of the group)
33+}
34+
35+interface Group {
36+ title: string;
37+ open: boolean;
38+ fields: FieldSpec[];
39+ selects?: SelectSpec[];
40+ note?: string | "grid" | "damping";
41+}
42+
43+const groups: Group[] = [
44+ {
45+ title: "Sheet",
46+ open: true,
47+ selects: [{ key: "shape", label: "Outline", options: [["rectangle", "Rectangle"], ["ellipse", "Ellipse (inscribed)"]] }],
48+ fields: [
49+ { key: "Lx", label: "Width", unit: "size", scale: 1, step: 0.5 },
50+ { key: "Ly", label: "Height", unit: "size", scale: 1, step: 0.5 },
51+ { key: "h", label: "Thickness", unit: "mm", scale: 1e-3, step: 0.01 },
52+ { key: "E", label: "Young's modulus", unit: "GPa", scale: 1e9, step: 1 },
53+ { key: "nu", label: "Poisson ratio", unit: "", scale: 1, step: 0.01 },
54+ { key: "rho", label: "Density", unit: "kg/m³", scale: 1, step: 10 },
55+ { key: "holeInset", label: "Hole inset", unit: "mm", scale: 1e-3, step: 1 },
56+ ],
57+ note: "grid",
58+ },
59+ {
60+ title: "Initial shape",
61+ open: true,
62+ selects: [{ key: "bowShape", label: "Bow shape", options: [["width", "Curved across the width"], ["height", "Curved across the height"], ["dish", "Dish (curved both ways)"], ["none", "Flat"]] }],
63+ fields: [{ key: "bowAmplitude", label: "Bow (center sag)", unit: "mm", scale: 1e-3, step: 0.5 }],
64+ note: "Real sheet metal is never flat. A bow of a few millimeters makes the sheet a shallow shell with an asymmetric stiffness; a negative push amount pushes against the bow.",
65+ },
66+ {
67+ title: "Excitation: push and release",
68+ open: true,
69+ fields: [
70+ { key: "pushAmount", label: "Push amount", unit: "mm", scale: 1e-3, step: 0.5 },
71+ { key: "fingerRadius", label: "Finger radius", unit: "mm", scale: 1e-3, step: 1 },
72+ { key: "rampTime", label: "Push ramp", unit: "s", scale: 1, step: 0.01 },
73+ { key: "holdMax", label: "Max settle time", unit: "s", scale: 1, step: 0.1 },
74+ { key: "releaseTime", label: "Release time", unit: "ms", scale: 1e-3, step: 1 },
75+ { key: "settleSigma0", label: "Settle damping σ₀", unit: "1/s", scale: 1, step: 10 },
76+ ],
77+ },
78+ {
79+ title: "Excitation: drumstick strike",
80+ open: false,
81+ fields: [
82+ { key: "strikeFracX", label: "Strike x (fraction of width)", unit: "", scale: 1, step: 0.05 },
83+ { key: "strikeFracY", label: "Strike y (fraction of height)", unit: "", scale: 1, step: 0.05 },
84+ { key: "strikeVelocity", label: "Stick speed", unit: "m/s", scale: 1, step: 0.5 },
85+ { key: "stickMass", label: "Stick mass", unit: "g", scale: 1e-3, step: 5 },
86+ { key: "contactStiffness", label: "Contact stiffness K", unit: "N/mᵅ", scale: 1, step: 1e8 },
87+ { key: "contactExponent", label: "Contact exponent α", unit: "", scale: 1, step: 0.1 },
88+ { key: "contactRadius", label: "Tip radius", unit: "mm", scale: 1e-3, step: 1 },
89+ ],
90+ },
91+ {
92+ title: "Losses",
93+ open: true,
94+ fields: [
95+ { key: "sigma0", label: "σ₀ (frequency independent)", unit: "1/s", scale: 1, step: 0.1, live: true },
96+ { key: "sigma1", label: "σ₁ (∝ frequency)", unit: "m²/s", scale: 1, step: 1e-5, live: true },
97+ { key: "sigmaIn", label: "In-plane loss", unit: "1/s", scale: 1, step: 5, live: true },
98+ ],
99+ note: "damping",
100+ },
101+ {
102+ title: "Listener",
103+ open: true,
104+ fields: [
105+ { key: "listenerDist", label: "Distance", unit: "m", scale: 1, step: 0.1, live: true },
106+ { key: "listenerAzimuth", label: "Azimuth", unit: "deg", scale: 1, step: 5, live: true },
107+ { key: "listenerElevation", label: "Elevation", unit: "deg", scale: 1, step: 5, live: true },
108+ ],
109+ },
110+ {
111+ title: "Numerics",
112+ open: false,
113+ fields: [
114+ { key: "spacing", label: "Grid spacing", unit: "mm", scale: 1e-3, step: 0.5 },
115+ { key: "fs", label: "Audio sample rate", unit: "Hz", scale: 1, step: 100 },
116+ { key: "fps", label: "Movie frame rate", unit: "fps", scale: 1, step: 10 },
117+ { key: "settleKE", label: "Settled below KE", unit: "J", scale: 1, step: 1e-9 },
118+ ],
119+ },
120+];
121+
122+const SIZE_SCALE: Record<SizeUnit, number> = { in: 0.0254, mm: 1e-3 };
123+const enumKeys: (keyof PhysParams)[] = ["shape", "bowShape", "excitation"];
124+const enumValues: Record<string, string[]> = {
125+ shape: ["rectangle", "ellipse"],
126+ bowShape: ["none", "width", "height", "dish"],
127+ excitation: ["push", "strike"],
128+};
129+
130+export function formatNumber(v: number): string {
131+ if (v === 0) return "0";
132+ const a = Math.abs(v);
133+ if (a >= 1e5 || a < 1e-3) return v.toExponential(3);
134+ return String(+v.toPrecision(6));
135+}
136+
137+export class ParamForm {
138+ readonly params: PhysParams;
139+ sizeUnit: SizeUnit = "in";
140+ /** Called after any change; `live` tells whether the field may change during a run. */
141+ onChange: (key: keyof PhysParams, live: boolean) => void = () => {};
142+
143+ private inputs = new Map<keyof PhysParams, HTMLInputElement | HTMLSelectElement>();
144+ private sizeUnitSelect: HTMLSelectElement;
145+ private gridNote: HTMLElement | null = null;
146+ private dampingNote: HTMLElement | null = null;
147+ private liveKeys = new Set<keyof PhysParams>();
148+
149+ constructor(root: HTMLElement, params: PhysParams) {
150+ this.params = params;
151+ // Excitation selector at the top, outside the groups.
152+ this.sizeUnitSelect = document.createElement("select");
153+ root.appendChild(this.selectRow("Excitation", "excitation", [["push", "Push a corner and release"], ["strike", "Drumstick strike"]]));
154+ for (const g of groups) {
155+ const det = document.createElement("details");
156+ det.open = g.open;
157+ const sum = document.createElement("summary");
158+ sum.textContent = g.title;
159+ det.appendChild(sum);
160+ for (const s of g.selects ?? []) if (!s.before) det.appendChild(this.selectRow(s.label, s.key, s.options));
161+ for (const f of g.fields) {
162+ for (const s of g.selects ?? []) if (s.before === f.key) det.appendChild(this.selectRow(s.label, s.key, s.options));
163+ det.appendChild(this.fieldRow(f));
164+ if (f.live) this.liveKeys.add(f.key);
165+ }
166+ if (g.note) {
167+ const note = document.createElement("p");
168+ note.className = "sub";
169+ if (g.note === "grid") this.gridNote = note;
170+ else if (g.note === "damping") this.dampingNote = note;
171+ else note.textContent = g.note;
172+ det.appendChild(note);
173+ }
174+ root.appendChild(det);
175+ }
176+ this.updateNotes();
177+ }
178+
179+ private selectRow(label: string, key: keyof PhysParams, options: [string, string][]): HTMLElement {
180+ const row = document.createElement("div");
181+ row.className = "field";
182+ const l = document.createElement("label");
183+ l.textContent = label;
184+ const sel = document.createElement("select");
185+ for (const [v, t] of options) {
186+ const o = document.createElement("option");
187+ o.value = v; o.textContent = t;
188+ sel.appendChild(o);
189+ }
190+ sel.value = String(this.params[key]);
191+ sel.addEventListener("change", () => {
192+ (this.params as unknown as Record<string, string>)[key] = sel.value;
193+ this.updateNotes();
194+ this.onChange(key, false);
195+ });
196+ row.append(l, sel, document.createElement("span"));
197+ this.inputs.set(key, sel);
198+ return row;
199+ }
200+
201+ private fieldRow(f: FieldSpec): HTMLElement {
202+ const row = document.createElement("div");
203+ row.className = "field";
204+ const label = document.createElement("label");
205+ label.textContent = f.label;
206+ const inp = document.createElement("input");
207+ inp.type = "number";
208+ inp.step = String(f.step ?? 1);
209+ inp.value = formatNumber((this.params[f.key] as number) / this.scale(f));
210+ inp.addEventListener("change", () => {
211+ const v = parseFloat(inp.value);
212+ if (!isFinite(v)) return;
213+ (this.params as unknown as Record<string, number>)[f.key] = v * this.scale(f);
214+ this.updateNotes();
215+ this.onChange(f.key, !!f.live);
216+ });
217+ this.inputs.set(f.key, inp);
218+ const unit = document.createElement("span");
219+ unit.className = "unit";
220+ if (f.unit === "size") {
221+ // Width and height share one unit selector (inches by default), shown on the first of them.
222+ const sel = this.sizeUnitSelect;
223+ if (sel.options.length === 0) {
224+ for (const u of ["in", "mm"] as SizeUnit[]) {
225+ const o = document.createElement("option");
226+ o.value = u; o.textContent = u;
227+ sel.appendChild(o);
228+ }
229+ sel.value = this.sizeUnit;
230+ sel.addEventListener("change", () => {
231+ this.sizeUnit = sel.value as SizeUnit;
232+ for (const el of this.sizeUnitLabels) el.textContent = this.sizeUnit;
233+ this.refresh();
234+ });
235+ row.append(label, inp, sel);
236+ return row;
237+ }
238+ unit.textContent = this.sizeUnit;
239+ this.sizeUnitLabels.push(unit);
240+ } else {
241+ unit.textContent = f.unit;
242+ }
243+ row.append(label, inp, unit);
244+ return row;
245+ }
246+ private sizeUnitLabels: HTMLElement[] = [];
247+
248+ private scale(f: FieldSpec): number {
249+ return f.unit === "size" ? SIZE_SCALE[this.sizeUnit] : f.scale;
250+ }
251+
252+ /** Rewrites every input from the current parameter values. */
253+ refresh(): void {
254+ for (const g of groups) {
255+ for (const f of g.fields) {
256+ const inp = this.inputs.get(f.key) as HTMLInputElement;
257+ inp.value = formatNumber((this.params[f.key] as number) / this.scale(f));
258+ }
259+ }
260+ for (const k of enumKeys) (this.inputs.get(k) as HTMLSelectElement).value = String(this.params[k]);
261+ this.updateNotes();
262+ }
263+
264+ /** Enables or disables the fields that cannot change while a simulation runs. */
265+ setEnabled(enabled: boolean): void {
266+ for (const [key, el] of this.inputs) if (!this.liveKeys.has(key)) el.disabled = !enabled;
267+ }
268+
269+ /** Applies URL parameters: any PhysParams key in SI units. Returns the keys that were set. */
270+ applyUrl(q: URLSearchParams): (keyof PhysParams)[] {
271+ const set: (keyof PhysParams)[] = [];
272+ for (const [k, v] of q) {
273+ if (!(k in this.params)) continue;
274+ const key = k as keyof PhysParams;
275+ if (enumKeys.includes(key)) {
276+ if (enumValues[key].includes(v)) {
277+ (this.params as unknown as Record<string, string>)[key] = v;
278+ set.push(key);
279+ }
280+ } else {
281+ const num = parseFloat(v);
282+ if (isFinite(num)) {
283+ (this.params as unknown as Record<string, number>)[key] = num;
284+ set.push(key);
285+ }
286+ }
287+ }
288+ this.refresh();
289+ return set;
290+ }
291+
292+ gridInfo(): { grid: Grid; inside: number } {
293+ const grid = makeGrid(this.params);
294+ return { grid, inside: countInside(grid) };
295+ }
296+
297+ private updateNotes(): void {
298+ const p = this.params;
299+ if (this.gridNote) {
300+ const { grid, inside } = this.gridInfo();
301+ const m = makeMaterial(p);
302+ const f3 = dispersionFrequency(m, grid, 0.03), f10 = dispersionFrequency(m, grid, 0.1);
303+ const outline = p.shape === "ellipse" ? "ellipse with a stepped edge" : "rectangle";
304+ this.gridNote.textContent = `${outline}; ${grid.nx} × ${grid.ny} grid, ${inside} nodes in the sheet. Bending partials come out low by 3% near ${(f3 / 1e3).toFixed(1)} kHz and by 10% near ${(f10 / 1e3).toFixed(1)} kHz (grid dispersion).`;
305+ }
306+ if (this.dampingNote) {
307+ const m = makeMaterial(p);
308+ const f = (hz: number) => t60(m, p, 2 * Math.PI * hz).toFixed(2);
309+ this.dampingNote.textContent = `T60: ${f(100)} s at 100 Hz, ${f(1000)} s at 1 kHz, ${f(5000)} s at 5 kHz.`;
310+ }
311+ }
312+
313+ get excitation(): Excitation {
314+ return this.params.excitation;
315+ }
316+ get shape(): PlateShape {
317+ return this.params.shape;
318+ }
319+ get bowShape(): BowShape {
320+ return this.params.bowShape;
321+ }
322+}
tests/cpuRef.test.tsadded+227−0View file
@@ -0,0 +1,227 @@
1+// Tests of the spatial discretization on the float64 reference (no GPU needed).
2+
3+import { test } from "node:test";
4+import assert from "node:assert/strict";
5+import { accelerations, addDamping, fields, maxAngularFrequency, potentialEnergy, type Fields } from "../src/sim/cpuRef.ts";
6+import { type Grid, type Material, countInside, nodeWeight, withEllipseMask } from "../src/params.ts";
7+
8+// 0.79 mm steel with nu = 0.3, the value the reference tables use.
9+const mat: Material = { D: 8.66, C: 1.665e8, G: 5.91e7, nu: 0.3, rhoh: 6.241 };
10+
11+function rectangle(nx: number, ny: number, Lx = 0.3, Ly = 0.2): Grid {
12+ return { nx, ny, hx: Lx / (nx - 1), hy: Ly / (ny - 1), n: nx * ny };
13+}
14+
15+function lcg(seed: number): () => number {
16+ return () => ((seed = (seed * 1664525 + 1013904223) >>> 0) / 4294967296 - 0.5);
17+}
18+
19+function massOf(g: Grid, k: number): number {
20+ return mat.rhoh * g.hx * g.hy * nodeWeight(g, k % g.nx, Math.floor(k / g.nx));
21+}
22+
23+/** Eigenvalues of a symmetric matrix (row-major, n x n) by cyclic Jacobi. */
24+function symEig(A: Float64Array, n: number): number[] {
25+ const a = Float64Array.from(A);
26+ for (let sweep = 0; sweep < 60; sweep++) {
27+ let off = 0;
28+ for (let p = 0; p < n; p++) for (let q = p + 1; q < n; q++) off += a[p * n + q] ** 2;
29+ if (off < 1e-22) break;
30+ for (let p = 0; p < n; p++) {
31+ for (let q = p + 1; q < n; q++) {
32+ const apq = a[p * n + q];
33+ if (Math.abs(apq) < 1e-300) continue;
34+ const theta = (a[q * n + q] - a[p * n + p]) / (2 * apq);
35+ const t = Math.sign(theta || 1) / (Math.abs(theta) + Math.sqrt(theta * theta + 1));
36+ const c = 1 / Math.sqrt(t * t + 1), s = t * c;
37+ for (let k = 0; k < n; k++) {
38+ const akp = a[k * n + p], akq = a[k * n + q];
39+ a[k * n + p] = c * akp - s * akq;
40+ a[k * n + q] = s * akp + c * akq;
41+ }
42+ for (let k = 0; k < n; k++) {
43+ const apk = a[p * n + k], aqk = a[q * n + k];
44+ a[p * n + k] = c * apk - s * aqk;
45+ a[q * n + k] = s * apk + c * aqk;
46+ }
47+ }
48+ }
49+ }
50+ const ev: number[] = [];
51+ for (let k = 0; k < n; k++) ev.push(a[k * n + k]);
52+ return ev.sort((x, y) => x - y);
53+}
54+
55+/** Mass-symmetrized linear stiffness M^{-1/2} K M^{-1/2} of one block ("w" or "uv") over the nodes of the sheet. */
56+function stiffness(g: Grid, block: "w" | "uv"): { K: Float64Array; n: number } {
57+ const nodes: number[] = [];
58+ for (let k = 0; k < g.n; k++) if (massOf(g, k) > 0) nodes.push(k);
59+ const N = nodes.length;
60+ const n = block === "w" ? N : 2 * N;
61+ const K = new Float64Array(n * n);
62+ const x = fields(g.n), y = fields(g.n);
63+ const eps = 1e-9;
64+ for (let col = 0; col < n; col++) {
65+ x.u.fill(0); x.v.fill(0); x.w.fill(0);
66+ if (block === "w") x.w[nodes[col]] = eps;
67+ else if (col < N) x.u[nodes[col]] = eps;
68+ else x.v[nodes[col - N]] = eps;
69+ accelerations(g, mat, x, y);
70+ const mc = massOf(g, nodes[col % N]);
71+ for (let row = 0; row < n; row++) {
72+ const acc = block === "w" ? y.w[nodes[row]] : row < N ? y.u[nodes[row]] : y.v[nodes[row - N]];
73+ K[row * n + col] = (-acc / eps) * Math.sqrt(massOf(g, nodes[row % N]) / mc);
74+ }
75+ }
76+ return { K, n };
77+}
78+
79+function checkSymmetricWithThreeRigidModes(K: Float64Array, n: number): number[] {
80+ let asym = 0, scale = 0;
81+ for (let r = 0; r < n; r++) {
82+ for (let c = 0; c < n; c++) {
83+ asym = Math.max(asym, Math.abs(K[r * n + c] - K[c * n + r]));
84+ scale = Math.max(scale, Math.abs(K[r * n + c]));
85+ }
86+ }
87+ assert.ok(asym < 1e-6 * scale, `asymmetry ${asym} vs scale ${scale}`);
88+ const ev = symEig(K, n);
89+ assert.ok(Math.abs(ev[2]) < 1e-6 * ev[n - 1], "third rigid mode");
90+ assert.ok(ev[3] > 1e-4 * ev[n - 1], `fourth mode is not rigid (${ev[3]} vs ${ev[n - 1]})`);
91+ return ev;
92+}
93+
94+/** Largest relative mismatch between the forces and the centerd-difference gradient of the energy. */
95+function gradientMismatch(g: Grid, f: Fields, w0?: Float64Array): number {
96+ const acc = fields(g.n);
97+ accelerations(g, mat, f, acc, w0);
98+ const delta = 1e-7;
99+ let worst = 0;
100+ for (const [arr, accArr] of [[f.w, acc.w], [f.u, acc.u], [f.v, acc.v]] as const) {
101+ for (let k = 0; k < g.n; k += 5) {
102+ const mass = massOf(g, k);
103+ if (mass === 0) continue;
104+ const orig = arr[k];
105+ arr[k] = orig + delta;
106+ const Ep = potentialEnergy(g, mat, f, w0);
107+ arr[k] = orig - delta;
108+ const Em = potentialEnergy(g, mat, f, w0);
109+ arr[k] = orig;
110+ const Ffd = -(Ep - Em) / (2 * delta);
111+ const Fan = accArr[k] * mass;
112+ worst = Math.max(worst, Math.abs(Ffd - Fan) / (Math.abs(Fan) + 1e-3));
113+ }
114+ }
115+ return worst;
116+}
117+
118+/** A random state with w of about 25 thicknesses (strongly nonlinear) and a cylindrical 3 mm bow. */
119+function randomState(g: Grid, seed: number): { f: Fields; w0: Float64Array } {
120+ const rnd = lcg(seed);
121+ const f = fields(g.n);
122+ const w0 = new Float64Array(g.n);
123+ for (let k = 0; k < g.n; k++) {
124+ const x = (k % g.nx) / (g.nx - 1) - 0.5;
125+ w0[k] = 0.003 * (1 - 4 * x * x);
126+ f.w[k] = 0.02 * rnd();
127+ f.u[k] = 1e-4 * rnd();
128+ f.v[k] = 1e-4 * rnd();
129+ }
130+ return { f, w0 };
131+}
132+
133+test("bending stiffness is symmetric and has exactly three rigid modes", () => {
134+ const { K, n } = stiffness(rectangle(9, 7), "w");
135+ checkSymmetricWithThreeRigidModes(K, n);
136+});
137+
138+test("membrane stiffness is symmetric and has exactly three rigid modes", () => {
139+ const { K, n } = stiffness(rectangle(7, 6), "uv");
140+ checkSymmetricWithThreeRigidModes(K, n);
141+});
142+
143+test("free square plate frequencies converge to the reference values", () => {
144+ // Completely free square plate, nu = 0.3, lambda^2 = omega a^2 sqrt(rho h / D)
145+ // (Gorman 1978; Leissa's Ritz values are slightly higher). The scheme is
146+ // second order: the errors are about 0.8 percent at 20 cells per side and
147+ // 0.3 percent at 30 (the pair at 34.80 is degenerate).
148+ const ref = [13.468, 19.596, 24.271, 34.801, 34.801, 61.093, 61.093];
149+ const a = 0.3;
150+ const errors: number[][] = [];
151+ for (const N of [16, 21]) {
152+ const { K, n } = stiffness(rectangle(N, N, a, a), "w");
153+ const ev = symEig(K, n);
154+ const lam2 = ev.slice(3, 3 + ref.length).map((e) => Math.sqrt(e) * a * a * Math.sqrt(mat.rhoh / mat.D));
155+ errors.push(lam2.map((v, k) => Math.abs(v - ref[k]) / ref[k]));
156+ }
157+ for (let k = 0; k < ref.length; k++) {
158+ assert.ok(errors[1][k] < 0.03, `mode ${k}: error ${(errors[1][k] * 100).toFixed(2)}% at 20 cells`);
159+ assert.ok(errors[1][k] < 0.7 * errors[0][k], `mode ${k}: error does not shrink from 15 to 20 cells`);
160+ }
161+});
162+
163+test("forces are the negative gradient of the discrete energy, flat and bowed", () => {
164+ const g = rectangle(8, 7);
165+ const { f, w0 } = randomState(g, 7);
166+ assert.ok(gradientMismatch(g, f) < 1e-4, "flat");
167+ assert.ok(gradientMismatch(g, f, w0) < 1e-4, "bowed");
168+ // The stress-free shape itself carries no energy.
169+ assert.ok(potentialEnergy(g, mat, fields(g.n), w0) < 1e-12, "w0 is stress free");
170+});
171+
172+test("loss terms only remove energy", () => {
173+ const g = rectangle(9, 8);
174+ const rnd = lcg(3);
175+ const v = fields(g.n), acc = fields(g.n);
176+ for (let k = 0; k < g.n; k++) { v.u[k] = rnd(); v.v[k] = rnd(); v.w[k] = rnd(); }
177+ addDamping(g, 1.0, 1e-3, 5.0, v, acc);
178+ let power = 0;
179+ for (let k = 0; k < g.n; k++) power += massOf(g, k) * (v.u[k] * acc.u[k] + v.v[k] * acc.v[k] + v.w[k] * acc.w[k]);
180+ assert.ok(power < 0, `power ${power} should be negative`);
181+});
182+
183+test("maximum frequency estimate is close to the in-plane CFL estimate", () => {
184+ const g: Grid = { nx: 50, ny: 50, hx: 0.004, hy: 0.004, n: 2500 };
185+ const wmax = maxAngularFrequency(g, mat);
186+ const cp = Math.sqrt(mat.C / mat.rhoh);
187+ const cfl = (2 * cp * Math.SQRT2) / 0.004; // 2 c sqrt(1/hx^2 + 1/hy^2)
188+ assert.ok(wmax > 0.5 * cfl && wmax < 1.5 * cfl, `wmax ${wmax} vs cfl ${cfl}`);
189+});
190+
191+// ---------- Masked (elliptical) grids ----------
192+
193+test("elliptical sheet: both stiffness blocks are symmetric with exactly three rigid modes", () => {
194+ const g = withEllipseMask(rectangle(13, 10));
195+ assert.ok(countInside(g) < g.n && countInside(g) > g.n / 2, "mask removes some nodes");
196+ for (let k = 0; k < g.n; k++) {
197+ // Every node of the sheet carries mass (no whiskers left after the mask is trimmed).
198+ if (g.mask![k]) assert.ok(massOf(g, k) > 0, `node ${k} has no mass`);
199+ }
200+ const b = stiffness(g, "w");
201+ checkSymmetricWithThreeRigidModes(b.K, b.n);
202+ const m = stiffness(g, "uv");
203+ checkSymmetricWithThreeRigidModes(m.K, m.n);
204+});
205+
206+test("elliptical sheet: forces are the energy gradient, flat and bowed", () => {
207+ const g = withEllipseMask(rectangle(12, 9));
208+ const { f, w0 } = randomState(g, 5);
209+ assert.ok(gradientMismatch(g, f) < 1e-4, "flat");
210+ assert.ok(gradientMismatch(g, f, w0) < 1e-4, "bowed");
211+});
212+
213+test("free circular plate frequencies approach the reference values", () => {
214+ const nx = 25, a = 0.15; // radius; 24 cells across the diameter
215+ const g = withEllipseMask({ nx, ny: nx, hx: (2 * a) / (nx - 1), hy: (2 * a) / (nx - 1), n: nx * nx });
216+ const { K, n } = stiffness(g, "w");
217+ const ev = symEig(K, n);
218+ const lam2 = ev.slice(3, 12).map((e) => Math.sqrt(e) * a * a * Math.sqrt(mat.rhoh / mat.D));
219+ // Completely free circular plate, nu = 0.3 (Leissa, Table 2.4): modes (2,0) x2, (0,1),
220+ // (3,0) x2, (1,1) x2, (4,0) x2 by nodal diameters and circles. The staircase edge is
221+ // first-order accurate: about 2 percent here, under 1 percent at 60 cells across.
222+ const ref = [5.358, 5.358, 9.003, 12.439, 12.439, 20.475, 20.475, 21.835, 21.835];
223+ for (let k = 0; k < ref.length; k++) {
224+ const relErr = Math.abs(lam2[k] - ref[k]) / ref[k];
225+ assert.ok(relErr < 0.04, `mode ${k}: got ${lam2[k].toFixed(3)}, expected ${ref[k]} (err ${(relErr * 100).toFixed(2)}%)`);
226+ }
227+});
tests/gpu.test.tsadded+111−0View file
@@ -0,0 +1,111 @@
1+// GPU-versus-CPU test: serves the app with Vite's dev server, opens it in a
2+// headless Chromium with WebGPU through Playwright, and reads the result of the
3+// in-page self-test (src/selftest.ts), which compares one step of the WGSL
4+// kernel with the float64 reference.
5+
6+import { test } from "node:test";
7+import assert from "node:assert/strict";
8+import { fileURLToPath } from "node:url";
9+import { dirname, join } from "node:path";
10+import { createServer } from "vite";
11+import { chromium, type Page } from "playwright";
12+
13+const root = join(dirname(fileURLToPath(import.meta.url)), "..");
14+
15+interface SelfTestCase {
16+ name: string;
17+ accError: { u: number; v: number; w: number };
18+ positionError: number;
19+ nodes: number;
20+ dt: number;
21+}
22+interface SelfTestResult { done: boolean; passed?: boolean; cases?: SelfTestCase[]; error?: string }
23+
24+interface Status {
25+ phase: string;
26+ t: number;
27+ totalSteps: number;
28+ wmax: number;
29+ peakPressure: number;
30+ samples: number;
31+ frames: number;
32+ stepsPerSec: number;
33+ dt: number;
34+ nx: number;
35+ ny: number;
36+ stick: { maxForce: number } | null;
37+}
38+
39+const ARGS = ["--enable-unsafe-webgpu", "--enable-features=Vulkan", "--use-angle=vulkan", "--ignore-gpu-blocklist", "--no-sandbox"];
40+
41+/** Serves the app and opens `query` in a headless Chromium; `body` gets the page and the list of page errors. */
42+async function withPage(query: string, body: (page: Page, errors: string[]) => Promise<void>): Promise<void> {
43+ const server = await createServer({ root, logLevel: "silent", server: { port: 0, strictPort: false } });
44+ await server.listen();
45+ const url = server.resolvedUrls?.local[0];
46+ assert.ok(url, "dev server URL");
47+ const browser = await chromium.launch({ headless: true, args: ARGS });
48+ try {
49+ const page = await browser.newPage();
50+ const errors: string[] = [];
51+ page.on("pageerror", (e) => errors.push(e.message));
52+ page.on("console", (msg) => { if (msg.type() === "error" || msg.type() === "warning") errors.push(`console.${msg.type()}: ${msg.text()}`); });
53+ await page.goto(`${url}?${query}`);
54+ await body(page, errors);
55+ } finally {
56+ await browser.close();
57+ await server.close();
58+ }
59+}
60+
61+test("one GPU step matches the float64 reference on flat, bowed, and elliptical sheets", { timeout: 120_000 }, async () => {
62+ await withPage("selftest=1", async (page, errors) => {
63+ let result: SelfTestResult = { done: false };
64+ for (let i = 0; i < 120 && !result.done; i++) {
65+ await page.waitForTimeout(500);
66+ result = await page.evaluate<SelfTestResult>("window.steelsim ? window.steelsim.selftest() : { done: false }");
67+ }
68+ assert.deepEqual(errors, [], "page errors");
69+ assert.ok(result.done, "self-test did not finish");
70+ assert.equal(result.error, undefined, `self-test error: ${result.error}`);
71+ for (const c of result.cases ?? []) {
72+ console.log(` ${c.name}: ${c.nodes} nodes, dt ${(c.dt * 1e6).toFixed(2)} us, acc errors ${c.accError.u.toExponential(1)} / ${c.accError.v.toExponential(1)} / ${c.accError.w.toExponential(1)}, position ${c.positionError.toExponential(1)}`);
73+ }
74+ assert.ok(result.passed, `self-test failed: ${JSON.stringify(result.cases)}`);
75+ });
76+});
77+
78+/** Runs a short simulation to completion and returns its final status. */
79+async function runToCompletion(query: string): Promise<Status> {
80+ let st: Status | null = null;
81+ await withPage(`autostart=1&${query}`, async (page, errors) => {
82+ for (let i = 0; i < 600; i++) {
83+ await page.waitForTimeout(250);
84+ st = await page.evaluate<Status | null>("window.steelsim ? window.steelsim.status() : null");
85+ const noGpu = await page.evaluate<string | null>("window.steelsim.noGpu()");
86+ assert.equal(noGpu, null, `no WebGPU: ${noGpu}`);
87+ if (st && st.phase === "done") break;
88+ }
89+ assert.deepEqual(errors, [], "page errors");
90+ });
91+ assert.ok(st, "no status");
92+ return st!;
93+}
94+
95+test("a short push-and-release run completes with finite output", { timeout: 180_000 }, async () => {
96+ const st = await runToCompletion("spacing=0.012&duration=0.2&rampTime=0.05&holdMax=0.1");
97+ console.log(` push: ${st.nx}x${st.ny}, dt ${(st.dt * 1e6).toFixed(2)} us, ${(st.stepsPerSec / 1e3).toFixed(0)} k steps/s, peak ${st.peakPressure.toFixed(3)} Pa, max |w| ${(st.wmax * 1e3).toFixed(2)} mm`);
98+ assert.equal(st.phase, "done");
99+ assert.ok(st.samples >= 0.2 * 48000 - 1, "recorded the requested duration");
100+ assert.ok(isFinite(st.wmax) && st.wmax < 0.05, "displacement stays bounded");
101+ assert.ok(st.peakPressure > 1e-3 && isFinite(st.peakPressure), "the sheet radiates");
102+});
103+
104+test("a short drumstick strike on an elliptical sheet completes with finite output", { timeout: 180_000 }, async () => {
105+ const st = await runToCompletion("spacing=0.012&duration=0.2&excitation=strike&shape=ellipse");
106+ console.log(` strike: ${st.nx}x${st.ny}, ${(st.stepsPerSec / 1e3).toFixed(0)} k steps/s, peak ${st.peakPressure.toFixed(3)} Pa, max contact force ${st.stick?.maxForce.toFixed(1)} N`);
107+ assert.equal(st.phase, "done");
108+ assert.ok(isFinite(st.wmax) && st.wmax < 0.05, "displacement stays bounded");
109+ assert.ok(st.peakPressure > 1e-3, "the sheet radiates");
110+ assert.ok(st.stick && st.stick.maxForce > 1, "the stick hit the sheet");
111+});
tsconfig.jsonadded+19−0View file
@@ -0,0 +1,19 @@
1+{
2+ "compilerOptions": {
3+ "target": "ES2022",
4+ "module": "ESNext",
5+ "moduleResolution": "bundler",
6+ "strict": true,
7+ "noUnusedLocals": true,
8+ "noUnusedParameters": true,
9+ "noFallthroughCasesInSwitch": true,
10+ "isolatedModules": true,
11+ "verbatimModuleSyntax": true,
12+ "skipLibCheck": true,
13+ "allowImportingTsExtensions": true,
14+ "noEmit": true,
15+ "lib": ["ES2022", "DOM", "DOM.Iterable"],
16+ "types": ["@webgpu/types", "vite/client", "node"]
17+ },
18+ "include": ["src", "tests"]
19+}
vite.config.tsadded+6−0View file
@@ -0,0 +1,6 @@
1+import { defineConfig } from "vite";
2+
3+export default defineConfig({
4+ base: "./",
5+ build: { target: "es2022" },
6+});