Nonlinear steel sheet gong simulator with WebGPU
Von Karman plate with in-plane inertia, energy-based free-edge bending
stencil, Q4 membrane elements, symplectic Euler stepping on the GPU,
rigid-body projection, Rayleigh-integral listener audio, 2D/3D views,
spectrogram, movie scrubbing, WAV download, and discretization tests.
26 changed files+3949−0
.gitignoreadded+2−0View file
@@ -0,0 +1,2 @@
1+node_modules/
2+dist/
README.mdadded+34−0View file
@@ -0,0 +1,34 @@
1+# steelsim21
2+
3+A browser simulation of the sound of a thin rectangular stainless steel sheet that hangs from a thread, is pushed by hand at the far corner, held until it settles, and released. 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+## Running
6+
7+Requires Node 20 or later and a browser with WebGPU (Chrome or Edge 113+; Firefox and Safari with WebGPU enabled).
8+
9+ npm install
10+ npm run dev # development server
11+ npm run build # static site in dist/
12+ npm test # discretization tests (Node, no GPU needed)
13+
14+## Using the simulator
15+
16+Press Start. The sheet is pushed at the bottom-right corner over the ramp time, held with heavy artificial damping until it settles, and then released; the recording of audio and movie frames starts at the release. The view shows the transverse displacement as a color map (2D) or a lit height field (3D, drag to orbit and 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.
17+
18+- Plate, excitation, and numerics parameters take effect at the next Start. Damping, listener position, and the radiation model can be changed while a simulation runs.
19+- After the simulation completes, drag the time slider to explore the movie; Play starts audio from the slider position and moves the movie and the spectrogram cursor with it. Live returns to the running view.
20+- To extend a finished simulation, raise the duration and press Continue; the state is kept on the GPU.
21+- Any parameter can be preset in the URL in SI units, for example `?spacing=0.006&duration=2&autostart=1`.
22+
23+The default grid spacing of 4 mm (about 8.9k nodes) resolves bending waves to about 20 kHz. With in-plane inertia the stable time step is set by the in-plane wave speed, about 0.65 microseconds, so four seconds of sound take about six million steps. This runs at roughly 35k steps per second on an Intel Iris Xe integrated GPU (about three minutes for four seconds) and considerably faster on a discrete GPU. A spacing of 6 mm is about three times faster and adequate below 10 kHz.
24+
25+## Layout
26+
27+- `src/params.ts`: physical and numerical parameters, derived quantities.
28+- `src/sim/cpuRef.ts`: float64 reference implementation of the spatial discretization; used for the time-step estimate and the tests.
29+- `src/sim/shaders/`: WGSL kernels (time step, rigid-body projection, listener sum, display extraction).
30+- `src/sim/simulation.ts`: WebGPU driver (buffers, batch encoding, phases, readback).
31+- `src/render/`: 2D and 3D WebGPU rendering of the displacement field.
32+- `src/audio/`: FFT, spectrogram, WAV encoder, playback.
33+- `src/main.ts`: user interface.
34+- `tests/`: symmetry, rigid modes, energy-gradient consistency, and free-plate frequencies against the Leissa reference values.
THEORY.mdadded+76−0View file
@@ -0,0 +1,76 @@
1+# Theory and numerical method
2+
3+This note describes the model behind the simulator: a thin stainless steel sheet, suspended by a thread, pushed at a corner, held, and released. The sheet rings like a gong, which means that its vibration amplitude (about ten thicknesses here) is far beyond the range where a linear plate model is accurate. We therefore solve the geometrically nonlinear plate equations of Foppl and von Karman with the in-plane inertia retained, discretize them on a regular grid, and step them explicitly on the GPU.
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]$ and has thickness $h$, Young's modulus $E$, Poisson ratio $\nu$, and density $\rho$. Its motion is described by three displacement fields, the in-plane displacements $u(x,y,t)$ and $v(x,y,t)$ and the transverse displacement $w(x,y,t)$. With the membrane strains
8+
9+$$
10+\varepsilon_{xx} = u_x + \tfrac12 w_x^2, \qquad
11+\varepsilon_{yy} = v_y + \tfrac12 w_y^2, \qquad
12+\gamma_{xy} = u_y + v_x + w_x w_y,
13+$$
14+
15+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}$, where $C = Eh/(1-\nu^2)$ and $G = Eh/(2(1+\nu))$. The equations of motion are
16+
17+$$
18+\rho h\, u_{tt} = \partial_x N_{xx} + \partial_y N_{xy}, \qquad
19+\rho h\, v_{tt} = \partial_x N_{xy} + \partial_y N_{yy},
20+$$
21+
22+$$
23+\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 ,
24+$$
25+
26+with the bending stiffness $D = Eh^3/(12(1-\nu^2))$. The quadratic terms in the strains couple bending to stretching. They are what produce the gong-like behaviour: the pitch glide of a strongly struck plate, and the cascade of energy from the low modes that the push excites toward high frequencies. The in-plane inertia terms $\rho h u_{tt}$ and $\rho h v_{tt}$ are usually dropped in favour of an Airy stress function, on the grounds that in-plane waves are much faster than bending waves. We keep them, as requested, at the cost of a much smaller time step (Section 5).
27+
28+The loss terms follow the two-parameter model that is standard in plate sound synthesis. A mode of angular frequency $\omega$ decays at the rate $\sigma_0 + \sigma_1 k^2$ with $k^2 = \omega/\kappa$ and $\kappa = \sqrt{D/\rho h}$, so $\sigma_0$ sets the low-frequency decay and $\sigma_1$ makes high frequencies decay faster. The defaults ($\sigma_0 = 1.15\ \mathrm{s^{-1}}$, $\sigma_1 = 1.08\times10^{-4}\ \mathrm{m^2/s}$) give a 60 dB decay time of about 6 s at low frequencies and about 1 s near 10 kHz. These are not measured values; both parameters are exposed in the interface.
29+
30+## 2. Boundary conditions and suspension
31+
32+All four edges are 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 the natural boundary conditions of that energy. This avoids the intricate ghost-point formulas that free edges otherwise require, and it guarantees a symmetric stiffness operator.
33+
34+The thread through the hole near one corner is not modelled as a constraint. At audio frequencies the thread is a very soft spring compared with the local stiffness of the plate, so it does not pin the hole; it only fixes the plate's average position. 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, applied 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}$, which is negligible.
35+
36+The same projection defines what "holding" the sheet means during the push. The finger prescribes the displacement of a small patch at the far corner. The reaction to the finger force is distributed over the plate as an inertial load (this is the inertia-relief equilibrium of a free body), so the push produces pure bending of the prescribed amount. This is a modelling 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 about 1.5 mm.
37+
38+## 3. Excitation
39+
40+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 constraint is removed instantaneously, the damping returns to its physical value, the clock is set to zero, and the recording of audio and movie frames begins. An instantaneous release is a step in the finger force and therefore excites all modes at once, in the same way that a plucked string is excited.
41+
42+## 4. Spatial discretization
43+
44+The plate is covered by a regular grid of $n_x \times n_y$ nodes with spacings $h_x$ and $h_y$ (about 4 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 and its linear part has a symmetric stiffness matrix.
45+
46+**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 they exist (not on the edges perpendicular to the direction of differencing), and $w_{xy}$ by the four-point difference at cell centres, with trapezoidal quadrature weights at edges and corners. 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, which yields the free-edge conditions. The three rigid modes $w = 1, x, y$ lie exactly in the null space.
47+
48+**Membrane.** The cells are treated as bilinear (Q4) finite elements with 2 by 2 Gauss integration, using the von Karman strains above. This is the standard element for this purpose; one-point integration would be cheaper but admits hourglass modes, which the nonlinear coupling would excite. The force on a node from each of its four cells is computed once per cell within a GPU workgroup and gathered by the nodes.
49+
50+**Mass.** The mass matrix is lumped with trapezoidal weights: a full cell mass in the interior, one half on an edge, one quarter at a corner, so that the total mass is exact.
51+
52+**Validation.** The tests in `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, and that the first five frequencies of a completely free square plate agree with the reference values of Leissa and Gorman ($\lambda^2 = 13.47, 19.60, 24.27, 34.80, 34.80$) to within 1 percent at 20 cells per side.
53+
54+## 5. Time integration
55+
56+The system is advanced with the symplectic Euler (Stormer-Verlet) 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. Its 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 4 mm spacing, $\Delta t \approx 0.65\ \mu$s, roughly six times smaller than bending alone would require, and 4 s of sound takes about 6 million steps. The step is then rounded down so that an integer number of steps (32 at 4 mm and 48 kHz) makes one audio sample.
57+
58+Single precision needs one precaution. With $w \approx 10$ mm 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).
59+
60+The state is stored twice and the kernel reads one copy and writes the other. Each step is one compute dispatch of 16 by 16 workgroups; a per-step uniform, addressed by a dynamic offset, carries the prescribed corner displacement, the damping, and the audio bookkeeping, so a whole movie frame of steps is encoded in one pass.
61+
62+## 6. Radiation to the listener
63+
64+The default listener model is the Rayleigh integral for a baffled plate,
65+
66+$$
67+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),
68+$$
69+
70+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 acceleration 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 common delay is discarded; only the spread of delays across the plate (about 5 samples for a listener 1 m in front) is kept. The plate is unbaffled in reality, so the sound from the two faces partly cancels at low frequencies. The "dipole" option is a heuristic for this: each element is weighted by the cosine of its angle to the listener, and the signal is high-pass filtered with a first-order filter at $c/(2\sqrt{L_x L_y})$, about 460 Hz for the default sheet. An accurate treatment would require solving the exterior acoustic problem, which is outside the scope of this simulator.
71+
72+The pressure is reported in pascals, and the peak sound pressure level is displayed. Playback and the downloaded WAV file are normalized to their peak.
73+
74+## 7. Limitations
75+
76+The von Karman model assumes moderate rotations, so it is not valid for displacements that are large compared with the plate size; the corner push of 10 mm on a 300 by 460 mm sheet is well within its range. The loss model is empirical and its parameters are not calibrated against a real sheet. The free-edge treatment near the edges is only first-order accurate. Air loading, which lowers the frequencies of a thin plate slightly and adds radiation damping, is not included. 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 gong 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 gong</h1>
13+ <p class="sub">Geometrically nonlinear plate (von Karman, in-plane inertia retained), explicit WebGPU time stepping.</p>
14+ <div id="controls"></div>
15+ <div id="params"></div>
16+ <div id="info"></div>
17+ <p class="sub">See THEORY.md in the repository for the model and the numerical method.</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+1187−0View file
@@ -0,0 +1,1187 @@
1+{
2+ "name": "steelsim21",
3+ "version": "0.1.0",
4+ "lockfileVersion": 3,
5+ "requires": true,
6+ "packages": {
7+ "": {
8+ "name": "steelsim21",
9+ "version": "0.1.0",
10+ "devDependencies": {
11+ "@webgpu/types": "^0.1.60",
12+ "typescript": "^5.6.0",
13+ "vite": "^6.0.0"
14+ }
15+ },
16+ "node_modules/@esbuild/aix-ppc64": {
17+ "version": "0.25.12",
18+ "resolved": "https://registry.npmjs.org/@esbuild/aix-ppc64/-/aix-ppc64-0.25.12.tgz",
19+ "integrity": "sha512-Hhmwd6CInZ3dwpuGTF8fJG6yoWmsToE+vYgD4nytZVxcu1ulHpUQRAB1UJ8+N1Am3Mz4+xOByoQoSZf4D+CpkA==",
20+ "cpu": [
21+ "ppc64"
22+ ],
23+ "dev": true,
24+ "license": "MIT",
25+ "optional": true,
26+ "os": [
27+ "aix"
28+ ],
29+ "engines": {
30+ "node": ">=18"
31+ }
32+ },
33+ "node_modules/@esbuild/android-arm": {
34+ "version": "0.25.12",
35+ "resolved": "https://registry.npmjs.org/@esbuild/android-arm/-/android-arm-0.25.12.tgz",
36+ "integrity": "sha512-VJ+sKvNA/GE7Ccacc9Cha7bpS8nyzVv0jdVgwNDaR4gDMC/2TTRc33Ip8qrNYUcpkOHUT5OZ0bUcNNVZQ9RLlg==",
37+ "cpu": [
38+ "arm"
39+ ],
40+ "dev": true,
41+ "license": "MIT",
42+ "optional": true,
43+ "os": [
44+ "android"
45+ ],
46+ "engines": {
47+ "node": ">=18"
48+ }
49+ },
50+ "node_modules/@esbuild/android-arm64": {
51+ "version": "0.25.12",
52+ "resolved": "https://registry.npmjs.org/@esbuild/android-arm64/-/android-arm64-0.25.12.tgz",
53+ "integrity": "sha512-6AAmLG7zwD1Z159jCKPvAxZd4y/VTO0VkprYy+3N2FtJ8+BQWFXU+OxARIwA46c5tdD9SsKGZ/1ocqBS/gAKHg==",
54+ "cpu": [
55+ "arm64"
56+ ],
57+ "dev": true,
58+ "license": "MIT",
59+ "optional": true,
60+ "os": [
61+ "android"
62+ ],
63+ "engines": {
64+ "node": ">=18"
65+ }
66+ },
67+ "node_modules/@esbuild/android-x64": {
68+ "version": "0.25.12",
69+ "resolved": "https://registry.npmjs.org/@esbuild/android-x64/-/android-x64-0.25.12.tgz",
70+ "integrity": "sha512-5jbb+2hhDHx5phYR2By8GTWEzn6I9UqR11Kwf22iKbNpYrsmRB18aX/9ivc5cabcUiAT/wM+YIZ6SG9QO6a8kg==",
71+ "cpu": [
72+ "x64"
73+ ],
74+ "dev": true,
75+ "license": "MIT",
76+ "optional": true,
77+ "os": [
78+ "android"
79+ ],
80+ "engines": {
81+ "node": ">=18"
82+ }
83+ },
84+ "node_modules/@esbuild/darwin-arm64": {
85+ "version": "0.25.12",
86+ "resolved": "https://registry.npmjs.org/@esbuild/darwin-arm64/-/darwin-arm64-0.25.12.tgz",
87+ "integrity": "sha512-N3zl+lxHCifgIlcMUP5016ESkeQjLj/959RxxNYIthIg+CQHInujFuXeWbWMgnTo4cp5XVHqFPmpyu9J65C1Yg==",
88+ "cpu": [
89+ "arm64"
90+ ],
91+ "dev": true,
92+ "license": "MIT",
93+ "optional": true,
94+ "os": [
95+ "darwin"
96+ ],
97+ "engines": {
98+ "node": ">=18"
99+ }
100+ },
101+ "node_modules/@esbuild/darwin-x64": {
102+ "version": "0.25.12",
103+ "resolved": "https://registry.npmjs.org/@esbuild/darwin-x64/-/darwin-x64-0.25.12.tgz",
104+ "integrity": "sha512-HQ9ka4Kx21qHXwtlTUVbKJOAnmG1ipXhdWTmNXiPzPfWKpXqASVcWdnf2bnL73wgjNrFXAa3yYvBSd9pzfEIpA==",
105+ "cpu": [
106+ "x64"
107+ ],
108+ "dev": true,
109+ "license": "MIT",
110+ "optional": true,
111+ "os": [
112+ "darwin"
113+ ],
114+ "engines": {
115+ "node": ">=18"
116+ }
117+ },
118+ "node_modules/@esbuild/freebsd-arm64": {
119+ "version": "0.25.12",
120+ "resolved": "https://registry.npmjs.org/@esbuild/freebsd-arm64/-/freebsd-arm64-0.25.12.tgz",
121+ "integrity": "sha512-gA0Bx759+7Jve03K1S0vkOu5Lg/85dou3EseOGUes8flVOGxbhDDh/iZaoek11Y8mtyKPGF3vP8XhnkDEAmzeg==",
122+ "cpu": [
123+ "arm64"
124+ ],
125+ "dev": true,
126+ "license": "MIT",
127+ "optional": true,
128+ "os": [
129+ "freebsd"
130+ ],
131+ "engines": {
132+ "node": ">=18"
133+ }
134+ },
135+ "node_modules/@esbuild/freebsd-x64": {
136+ "version": "0.25.12",
137+ "resolved": "https://registry.npmjs.org/@esbuild/freebsd-x64/-/freebsd-x64-0.25.12.tgz",
138+ "integrity": "sha512-TGbO26Yw2xsHzxtbVFGEXBFH0FRAP7gtcPE7P5yP7wGy7cXK2oO7RyOhL5NLiqTlBh47XhmIUXuGciXEqYFfBQ==",
139+ "cpu": [
140+ "x64"
141+ ],
142+ "dev": true,
143+ "license": "MIT",
144+ "optional": true,
145+ "os": [
146+ "freebsd"
147+ ],
148+ "engines": {
149+ "node": ">=18"
150+ }
151+ },
152+ "node_modules/@esbuild/linux-arm": {
153+ "version": "0.25.12",
154+ "resolved": "https://registry.npmjs.org/@esbuild/linux-arm/-/linux-arm-0.25.12.tgz",
155+ "integrity": "sha512-lPDGyC1JPDou8kGcywY0YILzWlhhnRjdof3UlcoqYmS9El818LLfJJc3PXXgZHrHCAKs/Z2SeZtDJr5MrkxtOw==",
156+ "cpu": [
157+ "arm"
158+ ],
159+ "dev": true,
160+ "license": "MIT",
161+ "optional": true,
162+ "os": [
163+ "linux"
164+ ],
165+ "engines": {
166+ "node": ">=18"
167+ }
168+ },
169+ "node_modules/@esbuild/linux-arm64": {
170+ "version": "0.25.12",
171+ "resolved": "https://registry.npmjs.org/@esbuild/linux-arm64/-/linux-arm64-0.25.12.tgz",
172+ "integrity": "sha512-8bwX7a8FghIgrupcxb4aUmYDLp8pX06rGh5HqDT7bB+8Rdells6mHvrFHHW2JAOPZUbnjUpKTLg6ECyzvas2AQ==",
173+ "cpu": [
174+ "arm64"
175+ ],
176+ "dev": true,
177+ "license": "MIT",
178+ "optional": true,
179+ "os": [
180+ "linux"
181+ ],
182+ "engines": {
183+ "node": ">=18"
184+ }
185+ },
186+ "node_modules/@esbuild/linux-ia32": {
187+ "version": "0.25.12",
188+ "resolved": "https://registry.npmjs.org/@esbuild/linux-ia32/-/linux-ia32-0.25.12.tgz",
189+ "integrity": "sha512-0y9KrdVnbMM2/vG8KfU0byhUN+EFCny9+8g202gYqSSVMonbsCfLjUO+rCci7pM0WBEtz+oK/PIwHkzxkyharA==",
190+ "cpu": [
191+ "ia32"
192+ ],
193+ "dev": true,
194+ "license": "MIT",
195+ "optional": true,
196+ "os": [
197+ "linux"
198+ ],
199+ "engines": {
200+ "node": ">=18"
201+ }
202+ },
203+ "node_modules/@esbuild/linux-loong64": {
204+ "version": "0.25.12",
205+ "resolved": "https://registry.npmjs.org/@esbuild/linux-loong64/-/linux-loong64-0.25.12.tgz",
206+ "integrity": "sha512-h///Lr5a9rib/v1GGqXVGzjL4TMvVTv+s1DPoxQdz7l/AYv6LDSxdIwzxkrPW438oUXiDtwM10o9PmwS/6Z0Ng==",
207+ "cpu": [
208+ "loong64"
209+ ],
210+ "dev": true,
211+ "license": "MIT",
212+ "optional": true,
213+ "os": [
214+ "linux"
215+ ],
216+ "engines": {
217+ "node": ">=18"
218+ }
219+ },
220+ "node_modules/@esbuild/linux-mips64el": {
221+ "version": "0.25.12",
222+ "resolved": "https://registry.npmjs.org/@esbuild/linux-mips64el/-/linux-mips64el-0.25.12.tgz",
223+ "integrity": "sha512-iyRrM1Pzy9GFMDLsXn1iHUm18nhKnNMWscjmp4+hpafcZjrr2WbT//d20xaGljXDBYHqRcl8HnxbX6uaA/eGVw==",
224+ "cpu": [
225+ "mips64el"
226+ ],
227+ "dev": true,
228+ "license": "MIT",
229+ "optional": true,
230+ "os": [
231+ "linux"
232+ ],
233+ "engines": {
234+ "node": ">=18"
235+ }
236+ },
237+ "node_modules/@esbuild/linux-ppc64": {
238+ "version": "0.25.12",
239+ "resolved": "https://registry.npmjs.org/@esbuild/linux-ppc64/-/linux-ppc64-0.25.12.tgz",
240+ "integrity": "sha512-9meM/lRXxMi5PSUqEXRCtVjEZBGwB7P/D4yT8UG/mwIdze2aV4Vo6U5gD3+RsoHXKkHCfSxZKzmDssVlRj1QQA==",
241+ "cpu": [
242+ "ppc64"
243+ ],
244+ "dev": true,
245+ "license": "MIT",
246+ "optional": true,
247+ "os": [
248+ "linux"
249+ ],
250+ "engines": {
251+ "node": ">=18"
252+ }
253+ },
254+ "node_modules/@esbuild/linux-riscv64": {
255+ "version": "0.25.12",
256+ "resolved": "https://registry.npmjs.org/@esbuild/linux-riscv64/-/linux-riscv64-0.25.12.tgz",
257+ "integrity": "sha512-Zr7KR4hgKUpWAwb1f3o5ygT04MzqVrGEGXGLnj15YQDJErYu/BGg+wmFlIDOdJp0PmB0lLvxFIOXZgFRrdjR0w==",
258+ "cpu": [
259+ "riscv64"
260+ ],
261+ "dev": true,
262+ "license": "MIT",
263+ "optional": true,
264+ "os": [
265+ "linux"
266+ ],
267+ "engines": {
268+ "node": ">=18"
269+ }
270+ },
271+ "node_modules/@esbuild/linux-s390x": {
272+ "version": "0.25.12",
273+ "resolved": "https://registry.npmjs.org/@esbuild/linux-s390x/-/linux-s390x-0.25.12.tgz",
274+ "integrity": "sha512-MsKncOcgTNvdtiISc/jZs/Zf8d0cl/t3gYWX8J9ubBnVOwlk65UIEEvgBORTiljloIWnBzLs4qhzPkJcitIzIg==",
275+ "cpu": [
276+ "s390x"
277+ ],
278+ "dev": true,
279+ "license": "MIT",
280+ "optional": true,
281+ "os": [
282+ "linux"
283+ ],
284+ "engines": {
285+ "node": ">=18"
286+ }
287+ },
288+ "node_modules/@esbuild/linux-x64": {
289+ "version": "0.25.12",
290+ "resolved": "https://registry.npmjs.org/@esbuild/linux-x64/-/linux-x64-0.25.12.tgz",
291+ "integrity": "sha512-uqZMTLr/zR/ed4jIGnwSLkaHmPjOjJvnm6TVVitAa08SLS9Z0VM8wIRx7gWbJB5/J54YuIMInDquWyYvQLZkgw==",
292+ "cpu": [
293+ "x64"
294+ ],
295+ "dev": true,
296+ "license": "MIT",
297+ "optional": true,
298+ "os": [
299+ "linux"
300+ ],
301+ "engines": {
302+ "node": ">=18"
303+ }
304+ },
305+ "node_modules/@esbuild/netbsd-arm64": {
306+ "version": "0.25.12",
307+ "resolved": "https://registry.npmjs.org/@esbuild/netbsd-arm64/-/netbsd-arm64-0.25.12.tgz",
308+ "integrity": "sha512-xXwcTq4GhRM7J9A8Gv5boanHhRa/Q9KLVmcyXHCTaM4wKfIpWkdXiMog/KsnxzJ0A1+nD+zoecuzqPmCRyBGjg==",
309+ "cpu": [
310+ "arm64"
311+ ],
312+ "dev": true,
313+ "license": "MIT",
314+ "optional": true,
315+ "os": [
316+ "netbsd"
317+ ],
318+ "engines": {
319+ "node": ">=18"
320+ }
321+ },
322+ "node_modules/@esbuild/netbsd-x64": {
323+ "version": "0.25.12",
324+ "resolved": "https://registry.npmjs.org/@esbuild/netbsd-x64/-/netbsd-x64-0.25.12.tgz",
325+ "integrity": "sha512-Ld5pTlzPy3YwGec4OuHh1aCVCRvOXdH8DgRjfDy/oumVovmuSzWfnSJg+VtakB9Cm0gxNO9BzWkj6mtO1FMXkQ==",
326+ "cpu": [
327+ "x64"
328+ ],
329+ "dev": true,
330+ "license": "MIT",
331+ "optional": true,
332+ "os": [
333+ "netbsd"
334+ ],
335+ "engines": {
336+ "node": ">=18"
337+ }
338+ },
339+ "node_modules/@esbuild/openbsd-arm64": {
340+ "version": "0.25.12",
341+ "resolved": "https://registry.npmjs.org/@esbuild/openbsd-arm64/-/openbsd-arm64-0.25.12.tgz",
342+ "integrity": "sha512-fF96T6KsBo/pkQI950FARU9apGNTSlZGsv1jZBAlcLL1MLjLNIWPBkj5NlSz8aAzYKg+eNqknrUJ24QBybeR5A==",
343+ "cpu": [
344+ "arm64"
345+ ],
346+ "dev": true,
347+ "license": "MIT",
348+ "optional": true,
349+ "os": [
350+ "openbsd"
351+ ],
352+ "engines": {
353+ "node": ">=18"
354+ }
355+ },
356+ "node_modules/@esbuild/openbsd-x64": {
357+ "version": "0.25.12",
358+ "resolved": "https://registry.npmjs.org/@esbuild/openbsd-x64/-/openbsd-x64-0.25.12.tgz",
359+ "integrity": "sha512-MZyXUkZHjQxUvzK7rN8DJ3SRmrVrke8ZyRusHlP+kuwqTcfWLyqMOE3sScPPyeIXN/mDJIfGXvcMqCgYKekoQw==",
360+ "cpu": [
361+ "x64"
362+ ],
363+ "dev": true,
364+ "license": "MIT",
365+ "optional": true,
366+ "os": [
367+ "openbsd"
368+ ],
369+ "engines": {
370+ "node": ">=18"
371+ }
372+ },
373+ "node_modules/@esbuild/openharmony-arm64": {
374+ "version": "0.25.12",
375+ "resolved": "https://registry.npmjs.org/@esbuild/openharmony-arm64/-/openharmony-arm64-0.25.12.tgz",
376+ "integrity": "sha512-rm0YWsqUSRrjncSXGA7Zv78Nbnw4XL6/dzr20cyrQf7ZmRcsovpcRBdhD43Nuk3y7XIoW2OxMVvwuRvk9XdASg==",
377+ "cpu": [
378+ "arm64"
379+ ],
380+ "dev": true,
381+ "license": "MIT",
382+ "optional": true,
383+ "os": [
384+ "openharmony"
385+ ],
386+ "engines": {
387+ "node": ">=18"
388+ }
389+ },
390+ "node_modules/@esbuild/sunos-x64": {
391+ "version": "0.25.12",
392+ "resolved": "https://registry.npmjs.org/@esbuild/sunos-x64/-/sunos-x64-0.25.12.tgz",
393+ "integrity": "sha512-3wGSCDyuTHQUzt0nV7bocDy72r2lI33QL3gkDNGkod22EsYl04sMf0qLb8luNKTOmgF/eDEDP5BFNwoBKH441w==",
394+ "cpu": [
395+ "x64"
396+ ],
397+ "dev": true,
398+ "license": "MIT",
399+ "optional": true,
400+ "os": [
401+ "sunos"
402+ ],
403+ "engines": {
404+ "node": ">=18"
405+ }
406+ },
407+ "node_modules/@esbuild/win32-arm64": {
408+ "version": "0.25.12",
409+ "resolved": "https://registry.npmjs.org/@esbuild/win32-arm64/-/win32-arm64-0.25.12.tgz",
410+ "integrity": "sha512-rMmLrur64A7+DKlnSuwqUdRKyd3UE7oPJZmnljqEptesKM8wx9J8gx5u0+9Pq0fQQW8vqeKebwNXdfOyP+8Bsg==",
411+ "cpu": [
412+ "arm64"
413+ ],
414+ "dev": true,
415+ "license": "MIT",
416+ "optional": true,
417+ "os": [
418+ "win32"
419+ ],
420+ "engines": {
421+ "node": ">=18"
422+ }
423+ },
424+ "node_modules/@esbuild/win32-ia32": {
425+ "version": "0.25.12",
426+ "resolved": "https://registry.npmjs.org/@esbuild/win32-ia32/-/win32-ia32-0.25.12.tgz",
427+ "integrity": "sha512-HkqnmmBoCbCwxUKKNPBixiWDGCpQGVsrQfJoVGYLPT41XWF8lHuE5N6WhVia2n4o5QK5M4tYr21827fNhi4byQ==",
428+ "cpu": [
429+ "ia32"
430+ ],
431+ "dev": true,
432+ "license": "MIT",
433+ "optional": true,
434+ "os": [
435+ "win32"
436+ ],
437+ "engines": {
438+ "node": ">=18"
439+ }
440+ },
441+ "node_modules/@esbuild/win32-x64": {
442+ "version": "0.25.12",
443+ "resolved": "https://registry.npmjs.org/@esbuild/win32-x64/-/win32-x64-0.25.12.tgz",
444+ "integrity": "sha512-alJC0uCZpTFrSL0CCDjcgleBXPnCrEAhTBILpeAp7M/OFgoqtAetfBzX0xM00MUsVVPpVjlPuMbREqnZCXaTnA==",
445+ "cpu": [
446+ "x64"
447+ ],
448+ "dev": true,
449+ "license": "MIT",
450+ "optional": true,
451+ "os": [
452+ "win32"
453+ ],
454+ "engines": {
455+ "node": ">=18"
456+ }
457+ },
458+ "node_modules/@napi-rs/lzma-linux-x64-gnu": {
459+ "version": "1.5.1",
460+ "resolved": "https://registry.npmjs.org/@napi-rs/lzma-linux-x64-gnu/-/lzma-linux-x64-gnu-1.5.1.tgz",
461+ "integrity": "sha512-oTXEIha4SsuXdTA4Iyskj0kpdx2yVXdhd75c2v3xGrHFfVMsbhTPZU/nMPL4sWKo4pBHm3aucLaqGlF696dTyQ==",
462+ "cpu": [
463+ "x64"
464+ ],
465+ "dev": true,
466+ "libc": [
467+ "glibc"
468+ ],
469+ "license": "MIT",
470+ "optional": true,
471+ "os": [
472+ "linux"
473+ ],
474+ "engines": {
475+ "node": "^22.20 || ^24.12 || >=25"
476+ }
477+ },
478+ "node_modules/@rollup/rollup-android-arm-eabi": {
479+ "version": "4.63.1",
480+ "resolved": "https://registry.npmjs.org/@rollup/rollup-android-arm-eabi/-/rollup-android-arm-eabi-4.63.1.tgz",
481+ "integrity": "sha512-UZ8sUxPTiHWYX9QNdJedb1kDZSpS1t/VPWBWGSgqHNi9w3Cu6IXvu2mzbhiTiPvtrqgTQJ+zqiAq2iPIPilpaQ==",
482+ "cpu": [
483+ "arm"
484+ ],
485+ "dev": true,
486+ "license": "MIT",
487+ "optional": true,
488+ "os": [
489+ "android"
490+ ]
491+ },
492+ "node_modules/@rollup/rollup-android-arm64": {
493+ "version": "4.63.1",
494+ "resolved": "https://registry.npmjs.org/@rollup/rollup-android-arm64/-/rollup-android-arm64-4.63.1.tgz",
495+ "integrity": "sha512-cQ4nFQABN5cDvDpbvJ7bMStCpnaVxynZrRMfUJYgxcIk9Sh54FIO1vtfkg0B69REjER77ioZ/ov+eAApx/KmLQ==",
496+ "cpu": [
497+ "arm64"
498+ ],
499+ "dev": true,
500+ "license": "MIT",
501+ "optional": true,
502+ "os": [
503+ "android"
504+ ]
505+ },
506+ "node_modules/@rollup/rollup-darwin-arm64": {
507+ "version": "4.63.1",
508+ "resolved": "https://registry.npmjs.org/@rollup/rollup-darwin-arm64/-/rollup-darwin-arm64-4.63.1.tgz",
509+ "integrity": "sha512-FQNqd1lRy/0QhDk3xeRIkSBiCpXCiDnZO3YLVdcDKN1UBiKToNftCzcXYNLshmPDUMlu2TdeS8tGcsU6f3YF1Q==",
510+ "cpu": [
511+ "arm64"
512+ ],
513+ "dev": true,
514+ "license": "MIT",
515+ "optional": true,
516+ "os": [
517+ "darwin"
518+ ]
519+ },
520+ "node_modules/@rollup/rollup-darwin-x64": {
521+ "version": "4.63.1",
522+ "resolved": "https://registry.npmjs.org/@rollup/rollup-darwin-x64/-/rollup-darwin-x64-4.63.1.tgz",
523+ "integrity": "sha512-pvD16V939D3CloK0+qikpGaxiPrDUXTe7Y5cWOMkMSy7m1cawa8EGy/kXYi/G/cKAC4HDAbSnzCIk1WmsoOKXg==",
524+ "cpu": [
525+ "x64"
526+ ],
527+ "dev": true,
528+ "license": "MIT",
529+ "optional": true,
530+ "os": [
531+ "darwin"
532+ ]
533+ },
534+ "node_modules/@rollup/rollup-freebsd-arm64": {
535+ "version": "4.63.1",
536+ "resolved": "https://registry.npmjs.org/@rollup/rollup-freebsd-arm64/-/rollup-freebsd-arm64-4.63.1.tgz",
537+ "integrity": "sha512-pcFGeL2345VwdTnJhA6zLbew+YgWB0qBG2+dMtXjCicf6+rm6kO6cOoh5VnTe0ZMrMRgRyuHmCJxZWrIdzYuOw==",
538+ "cpu": [
539+ "arm64"
540+ ],
541+ "dev": true,
542+ "license": "MIT",
543+ "optional": true,
544+ "os": [
545+ "freebsd"
546+ ]
547+ },
548+ "node_modules/@rollup/rollup-freebsd-x64": {
549+ "version": "4.63.1",
550+ "resolved": "https://registry.npmjs.org/@rollup/rollup-freebsd-x64/-/rollup-freebsd-x64-4.63.1.tgz",
551+ "integrity": "sha512-mRJlqSRulVzcKq/LKA6ICSIc3K/l4fzlVn/gePn2nXIHy8seRi5z/eeRE0d/XMBxcMldiXtQTSpRj0tkkC3g8Q==",
552+ "cpu": [
553+ "x64"
554+ ],
555+ "dev": true,
556+ "license": "MIT",
557+ "optional": true,
558+ "os": [
559+ "freebsd"
560+ ]
561+ },
562+ "node_modules/@rollup/rollup-linux-arm-gnueabihf": {
563+ "version": "4.63.1",
564+ "resolved": "https://registry.npmjs.org/@rollup/rollup-linux-arm-gnueabihf/-/rollup-linux-arm-gnueabihf-4.63.1.tgz",
565+ "integrity": "sha512-YDUNvVM85TI3g/1OpnqKP1h4NeW/j64DfWMf+G3M809xNk1bJSnpFp4sh83NpmVE5DXnkh8ULor4LTVZKoYLHw==",
566+ "cpu": [
567+ "arm"
568+ ],
569+ "dev": true,
570+ "libc": [
571+ "glibc"
572+ ],
573+ "license": "MIT",
574+ "optional": true,
575+ "os": [
576+ "linux"
577+ ]
578+ },
579+ "node_modules/@rollup/rollup-linux-arm-musleabihf": {
580+ "version": "4.63.1",
581+ "resolved": "https://registry.npmjs.org/@rollup/rollup-linux-arm-musleabihf/-/rollup-linux-arm-musleabihf-4.63.1.tgz",
582+ "integrity": "sha512-7Mcn71p9ZuQFAj+h+dhQXy/yeLePRS2yKRnmW1DijA9thKO5qap0GNOIQK4yQ6iP3SU0Mrb/yWo8h8vgRba8lw==",
583+ "cpu": [
584+ "arm"
585+ ],
586+ "dev": true,
587+ "libc": [
588+ "musl"
589+ ],
590+ "license": "MIT",
591+ "optional": true,
592+ "os": [
593+ "linux"
594+ ]
595+ },
596+ "node_modules/@rollup/rollup-linux-arm64-gnu": {
597+ "version": "4.63.1",
598+ "resolved": "https://registry.npmjs.org/@rollup/rollup-linux-arm64-gnu/-/rollup-linux-arm64-gnu-4.63.1.tgz",
599+ "integrity": "sha512-4YiLQTX6U4CSl0L9cluep9A9W6UmTfqBDc2/CH6wlu54pl4E7Jn3cOD8oxzvBDEGk/JMKgJ47C8g+radF7mwvg==",
600+ "cpu": [
601+ "arm64"
602+ ],
603+ "dev": true,
604+ "libc": [
605+ "glibc"
606+ ],
607+ "license": "MIT",
608+ "optional": true,
609+ "os": [
610+ "linux"
611+ ]
612+ },
613+ "node_modules/@rollup/rollup-linux-arm64-musl": {
614+ "version": "4.63.1",
615+ "resolved": "https://registry.npmjs.org/@rollup/rollup-linux-arm64-musl/-/rollup-linux-arm64-musl-4.63.1.tgz",
616+ "integrity": "sha512-2ra8F7w8OquwZN9z2/fKFnli69wa8PLwaVzRMIPGb13ByMJwC28Fbp8YcVGoUhlYMTt7j5j9bNgpysrN2UM+vw==",
617+ "cpu": [
618+ "arm64"
619+ ],
620+ "dev": true,
621+ "libc": [
622+ "musl"
623+ ],
624+ "license": "MIT",
625+ "optional": true,
626+ "os": [
627+ "linux"
628+ ]
629+ },
630+ "node_modules/@rollup/rollup-linux-loong64-gnu": {
631+ "version": "4.63.1",
632+ "resolved": "https://registry.npmjs.org/@rollup/rollup-linux-loong64-gnu/-/rollup-linux-loong64-gnu-4.63.1.tgz",
633+ "integrity": "sha512-Sy20ncyhjmBP0Ml+UvQbimjlk6VFgjW5uNP+qqwHB00mTE8Bl2C1TuHTlRwK2YoXeZbee5lP2XevBWVkAQAtSQ==",
634+ "cpu": [
635+ "loong64"
636+ ],
637+ "dev": true,
638+ "libc": [
639+ "glibc"
640+ ],
641+ "license": "MIT",
642+ "optional": true,
643+ "os": [
644+ "linux"
645+ ]
646+ },
647+ "node_modules/@rollup/rollup-linux-loong64-musl": {
648+ "version": "4.63.1",
649+ "resolved": "https://registry.npmjs.org/@rollup/rollup-linux-loong64-musl/-/rollup-linux-loong64-musl-4.63.1.tgz",
650+ "integrity": "sha512-noITLp8oNjYliPnGWmLyelIHwULGqbHloQHGw1rtxbWhTuWooRpnZarZQJ1y9EUC4szuCusCc+HEpUtxpIwYvA==",
651+ "cpu": [
652+ "loong64"
653+ ],
654+ "dev": true,
655+ "libc": [
656+ "musl"
657+ ],
658+ "license": "MIT",
659+ "optional": true,
660+ "os": [
661+ "linux"
662+ ]
663+ },
664+ "node_modules/@rollup/rollup-linux-ppc64-gnu": {
665+ "version": "4.63.1",
666+ "resolved": "https://registry.npmjs.org/@rollup/rollup-linux-ppc64-gnu/-/rollup-linux-ppc64-gnu-4.63.1.tgz",
667+ "integrity": "sha512-hlxxXd+F1mWiAcaFR7Sv9ZQT6m6UfI8+Vy/kFJzztq2pDMU/0wZ9sish0iszNZvsQDo8Gc0i5yuFEOz5dDf6fA==",
668+ "cpu": [
669+ "ppc64"
670+ ],
671+ "dev": true,
672+ "libc": [
673+ "glibc"
674+ ],
675+ "license": "MIT",
676+ "optional": true,
677+ "os": [
678+ "linux"
679+ ]
680+ },
681+ "node_modules/@rollup/rollup-linux-ppc64-musl": {
682+ "version": "4.63.1",
683+ "resolved": "https://registry.npmjs.org/@rollup/rollup-linux-ppc64-musl/-/rollup-linux-ppc64-musl-4.63.1.tgz",
684+ "integrity": "sha512-EF7OpqQTQ/BvGqLzUi4rEHuagCV9MugAUXSHemwPW5vxZ75RR+jxO/2j95Ph2dalMpFHSVECjRoioHZgA9zOYA==",
685+ "cpu": [
686+ "ppc64"
687+ ],
688+ "dev": true,
689+ "libc": [
690+ "musl"
691+ ],
692+ "license": "MIT",
693+ "optional": true,
694+ "os": [
695+ "linux"
696+ ]
697+ },
698+ "node_modules/@rollup/rollup-linux-riscv64-gnu": {
699+ "version": "4.63.1",
700+ "resolved": "https://registry.npmjs.org/@rollup/rollup-linux-riscv64-gnu/-/rollup-linux-riscv64-gnu-4.63.1.tgz",
701+ "integrity": "sha512-wQO3JesW9PRkwlabQ27y7sPfVOOTLRG73I4F2UYHG5PXun3J9U3y+b7ezVKSYbsvSKGQ1k1cq8Qlun4C9kLt3w==",
702+ "cpu": [
703+ "riscv64"
704+ ],
705+ "dev": true,
706+ "libc": [
707+ "glibc"
708+ ],
709+ "license": "MIT",
710+ "optional": true,
711+ "os": [
712+ "linux"
713+ ]
714+ },
715+ "node_modules/@rollup/rollup-linux-riscv64-musl": {
716+ "version": "4.63.1",
717+ "resolved": "https://registry.npmjs.org/@rollup/rollup-linux-riscv64-musl/-/rollup-linux-riscv64-musl-4.63.1.tgz",
718+ "integrity": "sha512-ouAGwhO6wHRXdnOVCOsB0tRFkA7nhNB2Nwax6oECXN0YiN8EYUTBAOudADOB1PI+yDL61TeNx/u7MVCzksNbkQ==",
719+ "cpu": [
720+ "riscv64"
721+ ],
722+ "dev": true,
723+ "libc": [
724+ "musl"
725+ ],
726+ "license": "MIT",
727+ "optional": true,
728+ "os": [
729+ "linux"
730+ ]
731+ },
732+ "node_modules/@rollup/rollup-linux-s390x-gnu": {
733+ "version": "4.63.1",
734+ "resolved": "https://registry.npmjs.org/@rollup/rollup-linux-s390x-gnu/-/rollup-linux-s390x-gnu-4.63.1.tgz",
735+ "integrity": "sha512-q2R38Sn+1J8RxhfJ+T54wSWmyKXWec+9jgDfqO2AtArEqHO5R2aeayp5H5OYLr5UYDVGsVaZPEFUooMhYCdz5A==",
736+ "cpu": [
737+ "s390x"
738+ ],
739+ "dev": true,
740+ "libc": [
741+ "glibc"
742+ ],
743+ "license": "MIT",
744+ "optional": true,
745+ "os": [
746+ "linux"
747+ ]
748+ },
749+ "node_modules/@rollup/rollup-linux-x64-gnu": {
750+ "version": "4.63.1",
751+ "resolved": "https://registry.npmjs.org/@rollup/rollup-linux-x64-gnu/-/rollup-linux-x64-gnu-4.63.1.tgz",
752+ "integrity": "sha512-gfI5T24WLLuFfSKw7Go/zDXjAAV0fny0swTaDv+WjK7vqcw4cRhFfdsyKL1n+ukI+ooBxn3bVQnyrn06WpI50w==",
753+ "cpu": [
754+ "x64"
755+ ],
756+ "dev": true,
757+ "libc": [
758+ "glibc"
759+ ],
760+ "license": "MIT",
761+ "optional": true,
762+ "os": [
763+ "linux"
764+ ]
765+ },
766+ "node_modules/@rollup/rollup-linux-x64-musl": {
767+ "version": "4.63.1",
768+ "resolved": "https://registry.npmjs.org/@rollup/rollup-linux-x64-musl/-/rollup-linux-x64-musl-4.63.1.tgz",
769+ "integrity": "sha512-4h6XqthmB4Hspji84wvgk+ElodTsGj+dbZqHJHHtKxj4mYq0ANSEEPX9ys3moJueqsRjwpaJYH7874Itwnj2ow==",
770+ "cpu": [
771+ "x64"
772+ ],
773+ "dev": true,
774+ "libc": [
775+ "musl"
776+ ],
777+ "license": "MIT",
778+ "optional": true,
779+ "os": [
780+ "linux"
781+ ]
782+ },
783+ "node_modules/@rollup/rollup-openbsd-x64": {
784+ "version": "4.63.1",
785+ "resolved": "https://registry.npmjs.org/@rollup/rollup-openbsd-x64/-/rollup-openbsd-x64-4.63.1.tgz",
786+ "integrity": "sha512-dlfCOa87o1VAYegLQ9EKilx2JCeRofiyPGhTCmqnuXZ6bMPiycO1rq1+sKoulAp7pGLIsTIw+1x5R+zgh5LhhA==",
787+ "cpu": [
788+ "x64"
789+ ],
790+ "dev": true,
791+ "license": "MIT",
792+ "optional": true,
793+ "os": [
794+ "openbsd"
795+ ]
796+ },
797+ "node_modules/@rollup/rollup-openharmony-arm64": {
798+ "version": "4.63.1",
799+ "resolved": "https://registry.npmjs.org/@rollup/rollup-openharmony-arm64/-/rollup-openharmony-arm64-4.63.1.tgz",
800+ "integrity": "sha512-cjkLbOlfcm3QGhMM1J5zaZjsw1GggbN6rw9UTSSRrPrR1KkcXnN7Uq9rPw34xImQ9VOY9GN+6u2Zj80B9ptkcw==",
801+ "cpu": [
802+ "arm64"
803+ ],
804+ "dev": true,
805+ "license": "MIT",
806+ "optional": true,
807+ "os": [
808+ "openharmony"
809+ ]
810+ },
811+ "node_modules/@rollup/rollup-win32-arm64-msvc": {
812+ "version": "4.63.1",
813+ "resolved": "https://registry.npmjs.org/@rollup/rollup-win32-arm64-msvc/-/rollup-win32-arm64-msvc-4.63.1.tgz",
814+ "integrity": "sha512-Li1KdUnWGE4N3e1F/B4RTB1ms+nG4WBgjByO46pkeBVX/2UBsY53xf5vK9WygVmnH3RwncIST7lkSdLSY6P9lg==",
815+ "cpu": [
816+ "arm64"
817+ ],
818+ "dev": true,
819+ "license": "MIT",
820+ "optional": true,
821+ "os": [
822+ "win32"
823+ ]
824+ },
825+ "node_modules/@rollup/rollup-win32-ia32-msvc": {
826+ "version": "4.63.1",
827+ "resolved": "https://registry.npmjs.org/@rollup/rollup-win32-ia32-msvc/-/rollup-win32-ia32-msvc-4.63.1.tgz",
828+ "integrity": "sha512-t4ZYOSoLTgwhuFMrmTMLx/+i1DQVK7HYqMc6kY46EApwi8X0nIVphzdNoThU3xt6n+N5urG1/gxBdCaKDLavfg==",
829+ "cpu": [
830+ "ia32"
831+ ],
832+ "dev": true,
833+ "license": "MIT",
834+ "optional": true,
835+ "os": [
836+ "win32"
837+ ]
838+ },
839+ "node_modules/@rollup/rollup-win32-x64-gnu": {
840+ "version": "4.63.1",
841+ "resolved": "https://registry.npmjs.org/@rollup/rollup-win32-x64-gnu/-/rollup-win32-x64-gnu-4.63.1.tgz",
842+ "integrity": "sha512-RgroPfMmKlD1RzSDxvwgcPiy2HNQKoYV7OmwIXDsk73uKW5t6B/V8KIy27SMv/FNXFo/oSBtWc9J0X7t91ezZg==",
843+ "cpu": [
844+ "x64"
845+ ],
846+ "dev": true,
847+ "license": "MIT",
848+ "optional": true,
849+ "os": [
850+ "win32"
851+ ]
852+ },
853+ "node_modules/@rollup/rollup-win32-x64-msvc": {
854+ "version": "4.63.1",
855+ "resolved": "https://registry.npmjs.org/@rollup/rollup-win32-x64-msvc/-/rollup-win32-x64-msvc-4.63.1.tgz",
856+ "integrity": "sha512-at8QVep6S3h5Y6gSbdGU06bRY5WJkf6WUduM9YtvYMbYhB1MOFfUgc6kehitQXzOtMSaT70q7f9ydPhpqu821w==",
857+ "cpu": [
858+ "x64"
859+ ],
860+ "dev": true,
861+ "license": "MIT",
862+ "optional": true,
863+ "os": [
864+ "win32"
865+ ]
866+ },
867+ "node_modules/@types/estree": {
868+ "version": "1.0.9",
869+ "resolved": "https://registry.npmjs.org/@types/estree/-/estree-1.0.9.tgz",
870+ "integrity": "sha512-GhdPgy1el4/ImP05X05Uw4cw2/M93BCUmnEvWZNStlCzEKME4Fkk+YpoA5OiHNQmoS7Cafb8Xa3Pya8m1Qrzeg==",
871+ "dev": true,
872+ "license": "MIT"
873+ },
874+ "node_modules/@webgpu/types": {
875+ "version": "0.1.72",
876+ "resolved": "https://registry.npmjs.org/@webgpu/types/-/types-0.1.72.tgz",
877+ "integrity": "sha512-0cF7RFM2edNoiIS1ODJp0/Gzv4/xSXhwoR0YCza+OWpJWtn4wmo9DvK91aLlH9+uUnwIriP7ZiC3WitmyhuzBw==",
878+ "dev": true,
879+ "license": "BSD-3-Clause"
880+ },
881+ "node_modules/esbuild": {
882+ "version": "0.25.12",
883+ "resolved": "https://registry.npmjs.org/esbuild/-/esbuild-0.25.12.tgz",
884+ "integrity": "sha512-bbPBYYrtZbkt6Os6FiTLCTFxvq4tt3JKall1vRwshA3fdVztsLAatFaZobhkBC8/BrPetoa0oksYoKXoG4ryJg==",
885+ "dev": true,
886+ "hasInstallScript": true,
887+ "license": "MIT",
888+ "bin": {
889+ "esbuild": "bin/esbuild"
890+ },
891+ "engines": {
892+ "node": ">=18"
893+ },
894+ "optionalDependencies": {
895+ "@esbuild/aix-ppc64": "0.25.12",
896+ "@esbuild/android-arm": "0.25.12",
897+ "@esbuild/android-arm64": "0.25.12",
898+ "@esbuild/android-x64": "0.25.12",
899+ "@esbuild/darwin-arm64": "0.25.12",
900+ "@esbuild/darwin-x64": "0.25.12",
901+ "@esbuild/freebsd-arm64": "0.25.12",
902+ "@esbuild/freebsd-x64": "0.25.12",
903+ "@esbuild/linux-arm": "0.25.12",
904+ "@esbuild/linux-arm64": "0.25.12",
905+ "@esbuild/linux-ia32": "0.25.12",
906+ "@esbuild/linux-loong64": "0.25.12",
907+ "@esbuild/linux-mips64el": "0.25.12",
908+ "@esbuild/linux-ppc64": "0.25.12",
909+ "@esbuild/linux-riscv64": "0.25.12",
910+ "@esbuild/linux-s390x": "0.25.12",
911+ "@esbuild/linux-x64": "0.25.12",
912+ "@esbuild/netbsd-arm64": "0.25.12",
913+ "@esbuild/netbsd-x64": "0.25.12",
914+ "@esbuild/openbsd-arm64": "0.25.12",
915+ "@esbuild/openbsd-x64": "0.25.12",
916+ "@esbuild/openharmony-arm64": "0.25.12",
917+ "@esbuild/sunos-x64": "0.25.12",
918+ "@esbuild/win32-arm64": "0.25.12",
919+ "@esbuild/win32-ia32": "0.25.12",
920+ "@esbuild/win32-x64": "0.25.12"
921+ }
922+ },
923+ "node_modules/fdir": {
924+ "version": "6.5.0",
925+ "resolved": "https://registry.npmjs.org/fdir/-/fdir-6.5.0.tgz",
926+ "integrity": "sha512-tIbYtZbucOs0BRGqPJkshJUYdL+SDH7dVM8gjy+ERp3WAUjLEFJE+02kanyHtwjWOnwrKYBiwAmM0p4kLJAnXg==",
927+ "dev": true,
928+ "license": "MIT",
929+ "engines": {
930+ "node": ">=12.0.0"
931+ },
932+ "peerDependencies": {
933+ "picomatch": "^3 || ^4"
934+ },
935+ "peerDependenciesMeta": {
936+ "picomatch": {
937+ "optional": true
938+ }
939+ }
940+ },
941+ "node_modules/fsevents": {
942+ "version": "2.3.3",
943+ "resolved": "https://registry.npmjs.org/fsevents/-/fsevents-2.3.3.tgz",
944+ "integrity": "sha512-5xoDfX+fL7faATnagmWPpbFtwh/R77WmMMqqHGS65C3vvB0YHrgF+B1YmZ3441tMj5n63k0212XNoJwzlhffQw==",
945+ "dev": true,
946+ "hasInstallScript": true,
947+ "license": "MIT",
948+ "optional": true,
949+ "os": [
950+ "darwin"
951+ ],
952+ "engines": {
953+ "node": "^8.16.0 || ^10.6.0 || >=11.0.0"
954+ }
955+ },
956+ "node_modules/nanoid": {
957+ "version": "3.3.18",
958+ "resolved": "https://registry.npmjs.org/nanoid/-/nanoid-3.3.18.tgz",
959+ "integrity": "sha512-DTg4MJbGMWkfi6VZFdNt2/caMbQy4Ou+Op/hJQvGEWcnVfoA1QA+xzRKAzw9jD6+GVOOeYr/mIcuDSdug6F6+w==",
960+ "dev": true,
961+ "funding": [
962+ {
963+ "type": "github",
964+ "url": "https://github.com/sponsors/ai"
965+ }
966+ ],
967+ "license": "MIT",
968+ "bin": {
969+ "nanoid": "bin/nanoid.cjs"
970+ },
971+ "engines": {
972+ "node": "^10 || ^12 || ^13.7 || ^14 || >=15.0.1"
973+ }
974+ },
975+ "node_modules/picocolors": {
976+ "version": "1.1.1",
977+ "resolved": "https://registry.npmjs.org/picocolors/-/picocolors-1.1.1.tgz",
978+ "integrity": "sha512-xceH2snhtb5M9liqDsmEw56le376mTZkEX/jEb/RxNFyegNul7eNslCXP9FDj/Lcu0X8KEyMceP2ntpaHrDEVA==",
979+ "dev": true,
980+ "license": "ISC"
981+ },
982+ "node_modules/picomatch": {
983+ "version": "4.0.7",
984+ "resolved": "https://registry.npmjs.org/picomatch/-/picomatch-4.0.7.tgz",
985+ "integrity": "sha512-qcJu88Q2IWqJsDD529JKMdwGm/dvInW4HvQnRwiH9JtihJvzGOscDtHE3x1pBKeUOTysQ8kVmLnJ2kJu7yhcGA==",
986+ "dev": true,
987+ "license": "MIT",
988+ "engines": {
989+ "node": ">=12"
990+ },
991+ "funding": {
992+ "url": "https://github.com/sponsors/jonschlinkert"
993+ }
994+ },
995+ "node_modules/postcss": {
996+ "version": "8.5.28",
997+ "resolved": "https://registry.npmjs.org/postcss/-/postcss-8.5.28.tgz",
998+ "integrity": "sha512-RRuzqDtt5Y9h3quz5hWhK+TPnsmVs6WwSU6LkJMeY4HstUEDuYTG8UJSdawMRzmzAtV+KEoG8N3Qg2qLy5vM/A==",
999+ "dev": true,
1000+ "funding": [
1001+ {
1002+ "type": "opencollective",
1003+ "url": "https://opencollective.com/postcss/"
1004+ },
1005+ {
1006+ "type": "tidelift",
1007+ "url": "https://tidelift.com/funding/github/npm/postcss"
1008+ },
1009+ {
1010+ "type": "github",
1011+ "url": "https://github.com/sponsors/ai"
1012+ }
1013+ ],
1014+ "license": "MIT",
1015+ "dependencies": {
1016+ "nanoid": "^3.3.18",
1017+ "picocolors": "^1.1.1",
1018+ "source-map-js": "^1.2.1"
1019+ },
1020+ "engines": {
1021+ "node": "^10 || ^12 || >=14"
1022+ }
1023+ },
1024+ "node_modules/rollup": {
1025+ "version": "4.63.1",
1026+ "resolved": "https://registry.npmjs.org/rollup/-/rollup-4.63.1.tgz",
1027+ "integrity": "sha512-3Df9jsstwhccuEfmAMi9l8XUh/GOkVObmFTU7CCVBysEbcOZLl84jCtaAZMcPiMz2EGKsATzQcU+Xr3n/wU6cg==",
1028+ "dev": true,
1029+ "license": "MIT",
1030+ "dependencies": {
1031+ "@types/estree": "1.0.9"
1032+ },
1033+ "bin": {
1034+ "rollup": "dist/bin/rollup"
1035+ },
1036+ "engines": {
1037+ "node": ">=18.0.0",
1038+ "npm": ">=8.0.0"
1039+ },
1040+ "optionalDependencies": {
1041+ "@napi-rs/lzma-linux-x64-gnu": "1.5.1",
1042+ "@rollup/rollup-android-arm-eabi": "4.63.1",
1043+ "@rollup/rollup-android-arm64": "4.63.1",
1044+ "@rollup/rollup-darwin-arm64": "4.63.1",
1045+ "@rollup/rollup-darwin-x64": "4.63.1",
1046+ "@rollup/rollup-freebsd-arm64": "4.63.1",
1047+ "@rollup/rollup-freebsd-x64": "4.63.1",
1048+ "@rollup/rollup-linux-arm-gnueabihf": "4.63.1",
1049+ "@rollup/rollup-linux-arm-musleabihf": "4.63.1",
1050+ "@rollup/rollup-linux-arm64-gnu": "4.63.1",
1051+ "@rollup/rollup-linux-arm64-musl": "4.63.1",
1052+ "@rollup/rollup-linux-loong64-gnu": "4.63.1",
1053+ "@rollup/rollup-linux-loong64-musl": "4.63.1",
1054+ "@rollup/rollup-linux-ppc64-gnu": "4.63.1",
1055+ "@rollup/rollup-linux-ppc64-musl": "4.63.1",
1056+ "@rollup/rollup-linux-riscv64-gnu": "4.63.1",
1057+ "@rollup/rollup-linux-riscv64-musl": "4.63.1",
1058+ "@rollup/rollup-linux-s390x-gnu": "4.63.1",
1059+ "@rollup/rollup-linux-x64-gnu": "4.63.1",
1060+ "@rollup/rollup-linux-x64-musl": "4.63.1",
1061+ "@rollup/rollup-openbsd-x64": "4.63.1",
1062+ "@rollup/rollup-openharmony-arm64": "4.63.1",
1063+ "@rollup/rollup-win32-arm64-msvc": "4.63.1",
1064+ "@rollup/rollup-win32-ia32-msvc": "4.63.1",
1065+ "@rollup/rollup-win32-x64-gnu": "4.63.1",
1066+ "@rollup/rollup-win32-x64-msvc": "4.63.1",
1067+ "fsevents": "~2.3.2"
1068+ }
1069+ },
1070+ "node_modules/source-map-js": {
1071+ "version": "1.2.1",
1072+ "resolved": "https://registry.npmjs.org/source-map-js/-/source-map-js-1.2.1.tgz",
1073+ "integrity": "sha512-UXWMKhLOwVKb728IUtQPXxfYU+usdybtUrK/8uGE8CQMvrhOpwvzDBwj0QhSL7MQc7vIsISBG8VQ8+IDQxpfQA==",
1074+ "dev": true,
1075+ "license": "BSD-3-Clause",
1076+ "engines": {
1077+ "node": ">=0.10.0"
1078+ }
1079+ },
1080+ "node_modules/tinyglobby": {
1081+ "version": "0.2.17",
1082+ "resolved": "https://registry.npmjs.org/tinyglobby/-/tinyglobby-0.2.17.tgz",
1083+ "integrity": "sha512-wXR/dYpcqKmfWpEdZjiKJOwCNFndD0DMnrW/cYjVGttEkBfVgcLFHoNrlj47mjOVic9yyNu65alsgF4NQyTa2g==",
1084+ "dev": true,
1085+ "license": "MIT",
1086+ "dependencies": {
1087+ "fdir": "^6.5.0",
1088+ "picomatch": "^4.0.4"
1089+ },
1090+ "engines": {
1091+ "node": ">=12.0.0"
1092+ },
1093+ "funding": {
1094+ "url": "https://github.com/sponsors/SuperchupuDev"
1095+ }
1096+ },
1097+ "node_modules/typescript": {
1098+ "version": "5.9.3",
1099+ "resolved": "https://registry.npmjs.org/typescript/-/typescript-5.9.3.tgz",
1100+ "integrity": "sha512-jl1vZzPDinLr9eUt3J/t7V6FgNEw9QjvBPdysz9KfQDD41fQrC2Y4vKQdiaUpFT4bXlb1RHhLpp8wtm6M5TgSw==",
1101+ "dev": true,
1102+ "license": "Apache-2.0",
1103+ "bin": {
1104+ "tsc": "bin/tsc",
1105+ "tsserver": "bin/tsserver"
1106+ },
1107+ "engines": {
1108+ "node": ">=14.17"
1109+ }
1110+ },
1111+ "node_modules/vite": {
1112+ "version": "6.4.3",
1113+ "resolved": "https://registry.npmjs.org/vite/-/vite-6.4.3.tgz",
1114+ "integrity": "sha512-NTKlcQjlAK7MlQoyb6LgaqHc8sso/pVyUJYWMws3jg21uTJw/LddqIFPcPqP6PzpgbIcZyKI85sFE4HBrQDA8A==",
1115+ "dev": true,
1116+ "license": "MIT",
1117+ "dependencies": {
1118+ "esbuild": "^0.25.0",
1119+ "fdir": "^6.4.4",
1120+ "picomatch": "^4.0.2",
1121+ "postcss": "^8.5.3",
1122+ "rollup": "^4.34.9",
1123+ "tinyglobby": "^0.2.13"
1124+ },
1125+ "bin": {
1126+ "vite": "bin/vite.js"
1127+ },
1128+ "engines": {
1129+ "node": "^18.0.0 || ^20.0.0 || >=22.0.0"
1130+ },
1131+ "funding": {
1132+ "url": "https://github.com/vitejs/vite?sponsor=1"
1133+ },
1134+ "optionalDependencies": {
1135+ "fsevents": "~2.3.3"
1136+ },
1137+ "peerDependencies": {
1138+ "@types/node": "^18.0.0 || ^20.0.0 || >=22.0.0",
1139+ "jiti": ">=1.21.0",
1140+ "less": "*",
1141+ "lightningcss": "^1.21.0",
1142+ "sass": "*",
1143+ "sass-embedded": "*",
1144+ "stylus": "*",
1145+ "sugarss": "*",
1146+ "terser": "^5.16.0",
1147+ "tsx": "^4.8.1",
1148+ "yaml": "^2.4.2"
1149+ },
1150+ "peerDependenciesMeta": {
1151+ "@types/node": {
1152+ "optional": true
1153+ },
1154+ "jiti": {
1155+ "optional": true
1156+ },
1157+ "less": {
1158+ "optional": true
1159+ },
1160+ "lightningcss": {
1161+ "optional": true
1162+ },
1163+ "sass": {
1164+ "optional": true
1165+ },
1166+ "sass-embedded": {
1167+ "optional": true
1168+ },
1169+ "stylus": {
1170+ "optional": true
1171+ },
1172+ "sugarss": {
1173+ "optional": true
1174+ },
1175+ "terser": {
1176+ "optional": true
1177+ },
1178+ "tsx": {
1179+ "optional": true
1180+ },
1181+ "yaml": {
1182+ "optional": true
1183+ }
1184+ }
1185+ }
1186+ }
1187+}
package.jsonadded+18−0View file
@@ -0,0 +1,18 @@
1+{
2+ "name": "steelsim21",
3+ "version": "0.1.0",
4+ "private": true,
5+ "type": "module",
6+ "description": "Browser-based nonlinear plate (Foppl-von Karman) sound simulation of a hanging stainless steel sheet, using WebGPU.",
7+ "scripts": {
8+ "dev": "vite",
9+ "build": "tsc --noEmit && vite build",
10+ "preview": "vite preview",
11+ "test": "node --test tests/*.test.ts"
12+ },
13+ "devDependencies": {
14+ "@webgpu/types": "^0.1.60",
15+ "typescript": "^5.6.0",
16+ "vite": "^6.0.0"
17+ }
18+}
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/player.tsadded+56−0View file
@@ -0,0 +1,56 @@
1+// Web Audio playback of a Float32Array signal with position reporting.
2+
3+export class Player {
4+ private ctx: AudioContext | null = null;
5+ private source: AudioBufferSourceNode | null = null;
6+ private startedAt = 0;
7+ private offset = 0;
8+ private length = 0;
9+ onEnded: (() => void) | null = null;
10+
11+ get playing(): boolean {
12+ return this.source !== null;
13+ }
14+
15+ /** Current playback position [s], or null when stopped. */
16+ position(): number | null {
17+ if (!this.source || !this.ctx) return null;
18+ return Math.min(this.length, this.offset + this.ctx.currentTime - this.startedAt);
19+ }
20+
21+ play(samples: Float32Array, fs: number, offset = 0, peak = 0.9): void {
22+ this.stop();
23+ if (samples.length === 0) return;
24+ if (!this.ctx) this.ctx = new AudioContext({ sampleRate: fs });
25+ const ctx = this.ctx;
26+ const buffer = ctx.createBuffer(1, samples.length, fs);
27+ const ch = buffer.getChannelData(0);
28+ let max = 0;
29+ for (let k = 0; k < samples.length; k++) max = Math.max(max, Math.abs(samples[k]));
30+ const gain = max > 0 ? peak / max : 0;
31+ for (let k = 0; k < samples.length; k++) ch[k] = samples[k] * gain;
32+ const src = ctx.createBufferSource();
33+ src.buffer = buffer;
34+ src.connect(ctx.destination);
35+ src.onended = () => {
36+ if (this.source === src) {
37+ this.source = null;
38+ this.onEnded?.();
39+ }
40+ };
41+ this.length = samples.length / fs;
42+ this.offset = Math.max(0, Math.min(offset, this.length));
43+ this.startedAt = ctx.currentTime;
44+ this.source = src;
45+ src.start(0, this.offset);
46+ void ctx.resume();
47+ }
48+
49+ stop(): void {
50+ if (this.source) {
51+ const s = this.source;
52+ this.source = null;
53+ try { s.stop(); } catch { /* already stopped */ }
54+ }
55+ }
56+}
src/audio/spectrogram.tsadded+186−0View file
@@ -0,0 +1,186 @@
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+ /** Draws the spectrogram; `duration` [s] fixes the time axis, `playhead` [s] draws a cursor. */
85+ draw(canvas: HTMLCanvasElement, fs: number, duration: number, playhead: number | null): void {
86+ const dpr = window.devicePixelRatio || 1;
87+ const W = Math.max(1, Math.floor(canvas.clientWidth * dpr));
88+ const H = Math.max(1, Math.floor(canvas.clientHeight * dpr));
89+ if (canvas.width !== W || canvas.height !== H) {
90+ canvas.width = W;
91+ canvas.height = H;
92+ }
93+ const ctx = canvas.getContext("2d")!;
94+ const left = 52 * dpr, bottom = 22 * dpr, top = 6 * dpr, right = 8 * dpr;
95+ const pw = Math.max(1, W - left - right), ph = Math.max(1, H - top - bottom);
96+ ctx.fillStyle = "#1c1e22";
97+ ctx.fillRect(0, 0, W, H);
98+ const { windowSize, hop, fmax, dbFloor, logFrequency } = this.opts;
99+ const binHz = fs / windowSize;
100+ const fmin = logFrequency ? 20 : 0;
101+ const img = ctx.createImageData(pw, ph);
102+ const data = img.data;
103+ const cols = this.columns.length;
104+ const colPerPx = (duration * fs) / hop / pw;
105+ const logMin = Math.log(fmin || 1), logMax = Math.log(fmax);
106+ for (let py = 0; py < ph; py++) {
107+ const frac = 1 - py / ph;
108+ const f = logFrequency ? Math.exp(logMin + frac * (logMax - logMin)) : fmin + frac * (fmax - fmin);
109+ const bin = Math.min(windowSize / 2 - 1, Math.max(0, Math.round(f / binHz)));
110+ for (let px = 0; px < pw; px++) {
111+ const c0 = Math.floor(px * colPerPx);
112+ const c1 = Math.max(c0 + 1, Math.floor((px + 1) * colPerPx));
113+ let v = -Infinity;
114+ for (let c = c0; c < c1 && c < cols; c++) v = Math.max(v, this.columns[c][bin]);
115+ const o = (py * pw + px) * 4;
116+ if (c0 >= cols) {
117+ data[o] = 28; data[o + 1] = 30; data[o + 2] = 34; data[o + 3] = 255;
118+ continue;
119+ }
120+ const t = this.peakDb > -200 ? Math.max(0, Math.min(1, (v - this.peakDb + dbFloor) / dbFloor)) : 0;
121+ const [r, g, b] = inferno(t);
122+ data[o] = r; data[o + 1] = g; data[o + 2] = b; data[o + 3] = 255;
123+ }
124+ }
125+ ctx.putImageData(img, left, top);
126+ // Axes.
127+ ctx.fillStyle = "#c8cbd0";
128+ ctx.font = `${11 * dpr}px system-ui, sans-serif`;
129+ ctx.textAlign = "right";
130+ ctx.textBaseline = "middle";
131+ const fticks = logFrequency ? [20, 50, 100, 200, 500, 1000, 2000, 5000, 10000, 20000].filter((f) => f >= fmin && f <= fmax) : niceTicks(0, fmax, 6);
132+ for (const f of fticks) {
133+ const frac = logFrequency ? (Math.log(f) - logMin) / (logMax - logMin) : (f - fmin) / (fmax - fmin);
134+ const y = top + (1 - frac) * ph;
135+ ctx.fillText(f >= 1000 ? `${f / 1000}k` : `${f}`, left - 4 * dpr, y);
136+ ctx.fillRect(left - 3 * dpr, y, 3 * dpr, 1);
137+ }
138+ ctx.textAlign = "center";
139+ ctx.textBaseline = "top";
140+ for (const t of niceTicks(0, duration, 8)) {
141+ const x = left + (t / duration) * pw;
142+ ctx.fillText(`${+t.toFixed(2)}`, x, top + ph + 4 * dpr);
143+ ctx.fillRect(x, top + ph, 1, 3 * dpr);
144+ }
145+ ctx.save();
146+ ctx.translate(12 * dpr, top + ph / 2);
147+ ctx.rotate(-Math.PI / 2);
148+ ctx.textBaseline = "middle";
149+ ctx.fillText("Hz", 0, 0);
150+ ctx.restore();
151+ if (playhead !== null && playhead >= 0 && playhead <= duration) {
152+ const x = left + (playhead / duration) * pw;
153+ ctx.fillStyle = "rgba(255,255,255,0.8)";
154+ ctx.fillRect(x, top, Math.max(1, dpr), ph);
155+ }
156+ }
157+}
158+
159+function niceTicks(lo: number, hi: number, count: number): number[] {
160+ const span = hi - lo;
161+ if (span <= 0) return [lo];
162+ const raw = span / count;
163+ const mag = 10 ** Math.floor(Math.log10(raw));
164+ const norm = raw / mag;
165+ const step = (norm < 1.5 ? 1 : norm < 3.5 ? 2 : norm < 7.5 ? 5 : 10) * mag;
166+ const ticks: number[] = [];
167+ for (let t = Math.ceil(lo / step) * step; t <= hi + 1e-9; t += step) ticks.push(t);
168+ return ticks;
169+}
170+
171+// Compact polynomial approximation of the inferno colormap.
172+function inferno(t: number): [number, number, number] {
173+ const c0 = [0.0002189403691192265, 0.001651004631001012, -0.01948089843709184];
174+ const c1 = [0.1065134194856116, 0.5639564367884091, 3.932712388889277];
175+ const c2 = [11.60249308247187, -3.972853965665698, -15.9423940062914];
176+ const c3 = [-41.70399613139459, 17.43639888205313, 44.35414519872813];
177+ const c4 = [77.162935699427, -33.40235894210092, -81.80730925738993];
178+ const c5 = [-71.31942824499214, 32.62606426397723, 73.20951985803202];
179+ const c6 = [25.13112622477341, -12.24266895238567, -23.07032500287172];
180+ const out: number[] = [];
181+ for (let k = 0; k < 3; k++) {
182+ const v = c0[k] + t * (c1[k] + t * (c2[k] + t * (c3[k] + t * (c4[k] + t * (c5[k] + t * c6[k])))));
183+ out.push(Math.max(0, Math.min(255, Math.round(v * 255))));
184+ }
185+ return [out[0], out[1], out[2]];
186+}
src/audio/wav.tsadded+28−0View file
@@ -0,0 +1,28 @@
1+// 16-bit PCM mono WAV encoder.
2+
3+export function encodeWav(samples: Float32Array, fs: number, peak = 0.99): Blob {
4+ let max = 0;
5+ for (let k = 0; k < samples.length; k++) max = Math.max(max, Math.abs(samples[k]));
6+ const gain = max > 0 ? peak / max : 0;
7+ const buf = new ArrayBuffer(44 + samples.length * 2);
8+ const v = new DataView(buf);
9+ const str = (o: number, s: string) => { for (let k = 0; k < s.length; k++) v.setUint8(o + k, s.charCodeAt(k)); };
10+ str(0, "RIFF");
11+ v.setUint32(4, 36 + samples.length * 2, true);
12+ str(8, "WAVE");
13+ str(12, "fmt ");
14+ v.setUint32(16, 16, true);
15+ v.setUint16(20, 1, true);
16+ v.setUint16(22, 1, true);
17+ v.setUint32(24, fs, true);
18+ v.setUint32(28, fs * 2, true);
19+ v.setUint16(32, 2, true);
20+ v.setUint16(34, 16, true);
21+ str(36, "data");
22+ v.setUint32(40, samples.length * 2, true);
23+ for (let k = 0; k < samples.length; k++) {
24+ const s = Math.max(-1, Math.min(1, samples[k] * gain));
25+ v.setInt16(44 + 2 * k, Math.round(s * 32767), true);
26+ }
27+ return new Blob([buf], { type: "audio/wav" });
28+}
src/main.tsadded+541−0View file
@@ -0,0 +1,541 @@
1+// UI wiring: parameters, run loop, live and movie views, playback, downloads.
2+
3+import { defaultParams, makeGrid, makeMaterial, t60, type PhysParams } from "./params.ts";
4+import { Simulation } from "./sim/simulation.ts";
5+import { Renderer, type ViewMode } from "./render/renderer.ts";
6+import { Spectrogram } from "./audio/spectrogram.ts";
7+import { Player } from "./audio/player.ts";
8+import { encodeWav } from "./audio/wav.ts";
9+
10+const $ = <T extends HTMLElement>(id: string) => document.getElementById(id) as T;
11+
12+// ---------- Parameter form ----------
13+
14+interface FieldSpec {
15+ key: keyof PhysParams;
16+ label: string;
17+ unit: string;
18+ scale: number; // UI value = SI value / scale
19+ step?: number;
20+ live?: boolean; // may be changed while a simulation runs
21+}
22+
23+const groups: { title: string; open: boolean; fields: FieldSpec[] }[] = [
24+ {
25+ title: "Plate",
26+ open: true,
27+ fields: [
28+ { key: "Lx", label: "Width", unit: "in", scale: 0.0254, step: 0.5 },
29+ { key: "Ly", label: "Height", unit: "in", scale: 0.0254, step: 0.5 },
30+ { key: "h", label: "Thickness", unit: "mm", scale: 1e-3, step: 0.01 },
31+ { key: "E", label: "Young's modulus", unit: "GPa", scale: 1e9, step: 1 },
32+ { key: "nu", label: "Poisson ratio", unit: "", scale: 1, step: 0.01 },
33+ { key: "rho", label: "Density", unit: "kg/m³", scale: 1, step: 10 },
34+ ],
35+ },
36+ {
37+ title: "Excitation",
38+ open: true,
39+ fields: [
40+ { key: "pushAmount", label: "Push at corner", unit: "mm", scale: 1e-3, step: 0.5 },
41+ { key: "fingerRadius", label: "Finger radius", unit: "mm", scale: 1e-3, step: 1 },
42+ { key: "rampTime", label: "Push ramp", unit: "s", scale: 1, step: 0.01 },
43+ { key: "holdMax", label: "Max settle time", unit: "s", scale: 1, step: 0.1 },
44+ { key: "settleSigma0", label: "Settle damping σ₀", unit: "1/s", scale: 1, step: 10 },
45+ { key: "holeInset", label: "Hole inset", unit: "mm", scale: 1e-3, step: 1 },
46+ ],
47+ },
48+ {
49+ title: "Damping",
50+ open: true,
51+ fields: [
52+ { key: "sigma0", label: "σ₀ (freq. independent)", unit: "1/s", scale: 1, step: 0.1, live: true },
53+ { key: "sigma1", label: "σ₁ (freq. dependent)", unit: "m²/s", scale: 1, step: 1e-5, live: true },
54+ ],
55+ },
56+ {
57+ title: "Listener",
58+ open: true,
59+ fields: [
60+ { key: "listenerDist", label: "Distance", unit: "m", scale: 1, step: 0.1, live: true },
61+ { key: "listenerAzimuth", label: "Azimuth", unit: "deg", scale: 1, step: 5, live: true },
62+ { key: "listenerElevation", label: "Elevation", unit: "deg", scale: 1, step: 5, live: true },
63+ ],
64+ },
65+ {
66+ title: "Numerics",
67+ open: false,
68+ fields: [
69+ { key: "spacing", label: "Grid spacing", unit: "mm", scale: 1e-3, step: 0.5 },
70+ { key: "fs", label: "Audio sample rate", unit: "Hz", scale: 1, step: 100 },
71+ { key: "fps", label: "Movie frame rate", unit: "fps", scale: 1, step: 10 },
72+ { key: "settleKE", label: "Settled below KE", unit: "J", scale: 1, step: 1e-9 },
73+ ],
74+ },
75+];
76+
77+const inputs = new Map<keyof PhysParams, HTMLInputElement | HTMLSelectElement>();
78+let params: PhysParams = { ...defaultParams };
79+
80+function buildForm(): void {
81+ const root = $("params");
82+ for (const g of groups) {
83+ const det = document.createElement("details");
84+ det.open = g.open;
85+ const sum = document.createElement("summary");
86+ sum.textContent = g.title;
87+ det.appendChild(sum);
88+ for (const f of g.fields) {
89+ const row = document.createElement("div");
90+ row.className = "field";
91+ const label = document.createElement("label");
92+ label.textContent = f.label;
93+ const inp = document.createElement("input");
94+ inp.type = "number";
95+ inp.step = String(f.step ?? 1);
96+ inp.value = fmt((params[f.key] as number) / f.scale);
97+ inp.addEventListener("change", () => onFieldChange(f, inp));
98+ const unit = document.createElement("span");
99+ unit.className = "unit";
100+ unit.textContent = f.unit;
101+ row.append(label, inp, unit);
102+ det.appendChild(row);
103+ inputs.set(f.key, inp);
104+ if (f.key === "listenerElevation") {
105+ const r2 = document.createElement("div");
106+ r2.className = "field";
107+ const l2 = document.createElement("label");
108+ l2.textContent = "Radiation model";
109+ const sel = document.createElement("select");
110+ for (const [v, t] of [["rayleigh", "Rayleigh (baffled)"], ["dipole", "Dipole (unbaffled, approx.)"]]) {
111+ const o = document.createElement("option");
112+ o.value = v; o.textContent = t;
113+ sel.appendChild(o);
114+ }
115+ sel.value = params.radiation;
116+ sel.addEventListener("change", () => {
117+ params.radiation = sel.value as PhysParams["radiation"];
118+ applyListener();
119+ });
120+ r2.append(l2, sel, document.createElement("span"));
121+ det.appendChild(r2);
122+ inputs.set("radiation", sel);
123+ }
124+ if (f.key === "sigma1") {
125+ const note = document.createElement("p");
126+ note.className = "sub";
127+ note.id = "dampNote";
128+ det.appendChild(note);
129+ }
130+ }
131+ root.appendChild(det);
132+ }
133+ updateDampingNote();
134+}
135+
136+function fmt(v: number): string {
137+ if (v === 0) return "0";
138+ const a = Math.abs(v);
139+ if (a >= 1e5 || a < 1e-3) return v.toExponential(3);
140+ return String(+v.toPrecision(6));
141+}
142+
143+function onFieldChange(f: FieldSpec, inp: HTMLInputElement): void {
144+ const v = parseFloat(inp.value);
145+ if (!isFinite(v)) return;
146+ (params as unknown as Record<string, number>)[f.key] = v * f.scale;
147+ if (f.key === "sigma0" || f.key === "sigma1") {
148+ updateDampingNote();
149+ sim?.setDamping(params.sigma0, params.sigma1);
150+ }
151+ if (f.key.startsWith("listener")) applyListener();
152+}
153+
154+function applyListener(): void {
155+ sim?.setListener(params.listenerDist, params.listenerAzimuth, params.listenerElevation, params.radiation);
156+}
157+
158+function updateDampingNote(): void {
159+ const m = makeMaterial(params);
160+ const f = (hz: number) => t60(m, params, 2 * Math.PI * hz).toFixed(2);
161+ $("dampNote").textContent = `T60: ${f(100)} s at 100 Hz, ${f(1000)} s at 1 kHz, ${f(5000)} s at 5 kHz, ${f(10000)} s at 10 kHz`;
162+}
163+
164+function setFormEnabled(enabled: boolean): void {
165+ for (const g of groups) for (const f of g.fields) if (!f.live) inputs.get(f.key)!.disabled = !enabled;
166+}
167+
168+// ---------- State ----------
169+
170+let device: GPUDevice | null = null;
171+let sim: Simulation | null = null;
172+let renderer: Renderer | null = null;
173+let running = false;
174+let paused = false;
175+let live = true;
176+let scrubFrame = 0;
177+let movieBuf: GPUBuffer | null = null;
178+let movieBufFrame = -1;
179+let spec = new Spectrogram({ windowSize: 2048, hop: 512, fmax: 10000, dbFloor: 80, logFrequency: false });
180+const player = new Player();
181+let lastStatsTime = 0, lastStatsSteps = 0, stepsPerSec = 0;
182+let duration = params.duration;
183+
184+// ---------- Controls ----------
185+
186+const btnStart = mkButton("Start", "primary");
187+const btnPause = mkButton("Pause");
188+const btnContinue = mkButton("Continue");
189+const durInput = document.createElement("input");
190+
191+function mkButton(text: string, cls = ""): HTMLButtonElement {
192+ const b = document.createElement("button");
193+ b.textContent = text;
194+ b.className = cls;
195+ return b;
196+}
197+
198+function buildControls(): void {
199+ const c = $("controls");
200+ const durRow = document.createElement("div");
201+ durRow.className = "field";
202+ durRow.style.width = "100%";
203+ const l = document.createElement("label");
204+ l.textContent = "Duration after release";
205+ durInput.type = "number";
206+ durInput.step = "0.5";
207+ durInput.min = "0.1";
208+ durInput.value = String(duration);
209+ durInput.addEventListener("change", () => {
210+ const v = parseFloat(durInput.value);
211+ if (isFinite(v) && v > 0) {
212+ duration = v;
213+ params.duration = v;
214+ sim?.setDuration(v);
215+ updateButtons();
216+ }
217+ });
218+ const u = document.createElement("span");
219+ u.className = "unit";
220+ u.textContent = "s";
221+ durRow.append(l, durInput, u);
222+ c.append(durRow, btnStart, btnPause, btnContinue);
223+ btnStart.addEventListener("click", () => void start());
224+ btnPause.addEventListener("click", () => {
225+ paused = !paused;
226+ btnPause.textContent = paused ? "Resume" : "Pause";
227+ if (!paused) void runLoop();
228+ });
229+ btnContinue.addEventListener("click", () => {
230+ if (!sim) return;
231+ sim.setDuration(duration);
232+ if (sim.phase !== "done") void runLoop();
233+ updateButtons();
234+ });
235+ updateButtons();
236+}
237+
238+function updateButtons(): void {
239+ const s = sim?.status();
240+ btnPause.disabled = !sim || !running || s?.phase === "done";
241+ btnContinue.disabled = !sim || running || !s || s.phase !== "done" || s.samples >= Math.round(duration * params.fs) - 1;
242+ if (sim && sim.phase === "done" && sim.samples < sim.targetSamples()) btnContinue.disabled = false;
243+}
244+
245+// ---------- View toolbars ----------
246+
247+const viewMode = { value: "2d" as ViewMode };
248+const scrub = document.createElement("input");
249+const timeLabel = document.createElement("span");
250+const btnLive = mkButton("Live");
251+const btnPlay = mkButton("Play");
252+const btnWav = mkButton("Download .wav");
253+let exagInput: HTMLInputElement;
254+let rangeInput: HTMLInputElement;
255+
256+function buildViewToolbar(): void {
257+ const tb = $("viewToolbar");
258+ const seg = document.createElement("span");
259+ seg.className = "seg";
260+ const b2 = mkButton("2D", "on"), b3 = mkButton("3D");
261+ b2.addEventListener("click", () => { viewMode.value = "2d"; b2.classList.add("on"); b3.classList.remove("on"); });
262+ b3.addEventListener("click", () => { viewMode.value = "3d"; b3.classList.add("on"); b2.classList.remove("on"); });
263+ seg.append(b2, b3);
264+ exagInput = numInput(4, 0.5);
265+ rangeInput = numInput(0, 0.5);
266+ const hint = document.createElement("span");
267+ hint.className = "sub";
268+ hint.style.margin = "0";
269+ hint.textContent = "3D: drag to orbit, wheel to zoom";
270+ tb.append(seg, labeled("Exaggeration ×", exagInput), labeled("Color range (mm, 0 = auto)", rangeInput), hint);
271+}
272+
273+function buildMovieToolbar(): void {
274+ const tb = $("movieToolbar");
275+ scrub.type = "range";
276+ scrub.min = "0"; scrub.max = "0"; scrub.value = "0";
277+ scrub.addEventListener("input", () => {
278+ live = false;
279+ scrubFrame = parseInt(scrub.value);
280+ if (player.playing) player.stop();
281+ updateLiveButton();
282+ });
283+ timeLabel.className = "time";
284+ btnLive.addEventListener("click", () => { live = true; updateLiveButton(); });
285+ btnPlay.addEventListener("click", () => {
286+ if (player.playing) { player.stop(); btnPlay.textContent = "Play"; return; }
287+ if (!sim || sim.samples === 0) return;
288+ const offset = live ? 0 : scrubFrame / sim.fps;
289+ live = false;
290+ player.play(sim.processedAudio(), sim.fs, offset);
291+ btnPlay.textContent = "Stop";
292+ updateLiveButton();
293+ });
294+ player.onEnded = () => { btnPlay.textContent = "Play"; };
295+ btnWav.addEventListener("click", () => {
296+ if (!sim || sim.samples === 0) return;
297+ const blob = encodeWav(sim.processedAudio(), sim.fs);
298+ const a = document.createElement("a");
299+ a.href = URL.createObjectURL(blob);
300+ a.download = `sheet_${sim.samples}samples_${sim.fs}Hz.wav`;
301+ a.click();
302+ setTimeout(() => URL.revokeObjectURL(a.href), 5000);
303+ });
304+ tb.append(btnPlay, btnLive, scrub, timeLabel, btnWav);
305+ updateLiveButton();
306+}
307+
308+function updateLiveButton(): void {
309+ btnLive.classList.toggle("on", live);
310+ btnLive.disabled = live;
311+}
312+
313+function buildSpecToolbar(): void {
314+ const tb = $("specToolbar");
315+ const win = document.createElement("select");
316+ for (const n of [512, 1024, 2048, 4096, 8192]) {
317+ const o = document.createElement("option");
318+ o.value = String(n); o.textContent = String(n);
319+ win.appendChild(o);
320+ }
321+ win.value = "2048";
322+ const hop = numInput(512, 64);
323+ const fmax = numInput(10000, 500);
324+ const floor = numInput(80, 5);
325+ const logf = document.createElement("input");
326+ logf.type = "checkbox";
327+ const apply = () => {
328+ if (!sim) return;
329+ spec.configure(
330+ { 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 },
331+ sim.processedAudio(), sim.samples,
332+ );
333+ };
334+ for (const el of [win, hop, fmax, floor, logf]) el.addEventListener("change", apply);
335+ tb.append(
336+ document.createTextNode("Spectrogram"),
337+ labeled("Window", win), labeled("Hop", hop), labeled("Max Hz", fmax), labeled("dB range", floor), labeled("Log freq.", logf),
338+ );
339+}
340+
341+function numInput(v: number, step: number): HTMLInputElement {
342+ const i = document.createElement("input");
343+ i.type = "number"; i.value = String(v); i.step = String(step);
344+ return i;
345+}
346+function labeled(text: string, el: HTMLElement): HTMLLabelElement {
347+ const l = document.createElement("label");
348+ l.append(text, el);
349+ return l;
350+}
351+
352+// ---------- Simulation control ----------
353+
354+async function start(): Promise<void> {
355+ if (!device) return;
356+ player.stop();
357+ running = false;
358+ await new Promise((r) => setTimeout(r, 0));
359+ sim?.destroy();
360+ spec.reset();
361+ live = true;
362+ paused = false;
363+ btnPause.textContent = "Pause";
364+ params.duration = duration;
365+ sim = new Simulation(device, params);
366+ renderer = new Renderer(device, $("view"), sim.grid, params);
367+ movieBuf?.destroy();
368+ movieBuf = device.createBuffer({ size: Math.max(16, sim.grid.n * 4), usage: GPUBufferUsage.STORAGE | GPUBufferUsage.COPY_DST });
369+ movieBufFrame = -1;
370+ scrub.max = "0"; scrub.value = "0";
371+ lastStatsTime = performance.now(); lastStatsSteps = 0; stepsPerSec = 0;
372+ setFormEnabled(false);
373+ updateLiveButton();
374+ void runLoop();
375+}
376+
377+async function runLoop(): Promise<void> {
378+ if (!sim || running) return;
379+ running = true;
380+ updateButtons();
381+ try {
382+ while (sim.phase !== "done" && !paused) {
383+ await sim.runBatch();
384+ const audio = sim.processedAudio();
385+ spec.update(audio, sim.samples);
386+ scrub.max = String(Math.max(0, sim.frames.length - 1));
387+ if (live) scrub.value = scrub.max;
388+ const now = performance.now();
389+ if (now - lastStatsTime > 500) {
390+ stepsPerSec = ((sim.totalSteps - lastStatsSteps) * 1000) / (now - lastStatsTime);
391+ lastStatsTime = now; lastStatsSteps = sim.totalSteps;
392+ }
393+ }
394+ } catch (e) {
395+ console.error(e);
396+ $("info").innerHTML = `<div class="err">Simulation error: ${String(e)}</div>`;
397+ }
398+ running = false;
399+ setFormEnabled(sim.phase === "done" && !paused ? true : false);
400+ if (sim.phase === "done") setFormEnabled(true);
401+ updateButtons();
402+}
403+
404+// ---------- Render loop ----------
405+
406+let frameErrorShown = false;
407+function frame(): void {
408+ requestAnimationFrame(frame);
409+ try {
410+ frameBody();
411+ } catch (e) {
412+ if (!frameErrorShown) { frameErrorShown = true; console.error("frame exception:", e); }
413+ }
414+}
415+
416+function frameBody(): void {
417+ if (!device || !sim || !renderer) return;
418+ const dev = device;
419+ renderer.mode = viewMode.value;
420+ renderer.exaggeration = parseFloat(exagInput.value) || 1;
421+ renderer.colorRange = (parseFloat(rangeInput.value) || 0) * 1e-3;
422+ const pos = player.position();
423+ if (pos !== null) {
424+ scrubFrame = Math.min(sim.frames.length - 1, Math.floor(pos * sim.fps));
425+ scrub.value = String(scrubFrame);
426+ }
427+ const enc = dev.createCommandEncoder();
428+ const st = sim.status();
429+ if (live || sim.frames.length === 0) {
430+ sim.encodeDisplay(enc);
431+ renderer.render(enc, sim.displayW, st.wmax);
432+ timeLabel.textContent = `${st.t >= 0 ? "" : "prep "}${st.t.toFixed(3)} s`;
433+ } else {
434+ const f = Math.max(0, Math.min(sim.frames.length - 1, scrubFrame));
435+ if (movieBufFrame !== f) {
436+ dev.queue.writeBuffer(movieBuf!, 0, sim.frames[f]);
437+ movieBufFrame = f;
438+ }
439+ let wmax = 0;
440+ const fr = sim.frames[f];
441+ for (let k = 0; k < fr.length; k += 7) wmax = Math.max(wmax, Math.abs(fr[k]));
442+ renderer.render(enc, movieBuf!, wmax);
443+ timeLabel.textContent = `${(f / sim.fps).toFixed(3)} s`;
444+ }
445+ dev.pushErrorScope("validation");
446+ dev.queue.submit([enc.finish()]);
447+ void dev.popErrorScope().then((e) => { if (e && !frameErrorShown) { frameErrorShown = true; console.error("frame validation error:", e.message); } });
448+ spec.draw($("spec"), sim.fs, Math.max(duration, sim.samples / sim.fs), pos ?? (live ? null : scrubFrame / sim.fps));
449+ updateInfo(st);
450+}
451+
452+function updateInfo(st: ReturnType<Simulation["status"]>): void {
453+ if (!sim) return;
454+ const g = sim.grid;
455+ const spl = st.peakPressure > 0 ? 20 * Math.log10(st.peakPressure / 2e-5) : -Infinity;
456+ const phase = { push: "pushing", hold: "settling (held)", recording: "ringing (recording)", done: "done" }[st.phase];
457+ const rows: [string, string][] = [
458+ ["Phase", phase],
459+ ["Time", `${st.t.toFixed(3)} s`],
460+ ["Grid", `${g.nx} × ${g.ny} nodes, ${(g.hx * 1e3).toFixed(2)} × ${(g.hy * 1e3).toFixed(2)} mm`],
461+ ["Time step", `${(sim.dt * 1e6).toFixed(3)} µs (${sim.stepsPerSample} per sample)`],
462+ ["Speed", `${(stepsPerSec / 1e3).toFixed(0)} k steps/s (${(stepsPerSec * sim.dt).toFixed(3)}× real time)`],
463+ ["Kinetic energy", `${st.ke.toExponential(2)} J`],
464+ ["Max |w|", `${(st.wmax * 1e3).toFixed(3)} mm`],
465+ ["Audio", `${(st.samples / sim.fs).toFixed(2)} s, ${st.frames} frames`],
466+ ["Peak pressure", st.peakPressure > 0 ? `${st.peakPressure.toFixed(3)} Pa (${spl.toFixed(1)} dB SPL)` : "–"],
467+ ];
468+ let html = rows.map(([a, b]) => `<div class="row"><span>${a}</span><span>${b}</span></div>`).join("");
469+ if (st.delayClamped) html += `<div class="warn">Listener is so close that some propagation delays were clamped.</div>`;
470+ if (!isFinite(st.wmax) || st.wmax > 1) html += `<div class="err">The solution has diverged. Reduce the time step by increasing the grid spacing, or reduce the push.</div>`;
471+ $("info").innerHTML = html;
472+}
473+
474+// ---------- Init ----------
475+
476+async function init(): Promise<void> {
477+ buildForm();
478+ buildControls();
479+ buildViewToolbar();
480+ buildMovieToolbar();
481+ buildSpecToolbar();
482+ if (!("gpu" in navigator)) {
483+ showNoGpu("WebGPU is not available in this browser. Use a recent Chrome or Edge (or Firefox / Safari with WebGPU enabled).");
484+ return;
485+ }
486+ const adapter = await navigator.gpu.requestAdapter({ powerPreference: "high-performance" });
487+ if (!adapter) {
488+ showNoGpu("No WebGPU adapter was found.");
489+ return;
490+ }
491+ device = await adapter.requestDevice();
492+ device.lost.then((info) => showNoGpu(`WebGPU device lost: ${info.message}`));
493+ device.addEventListener("uncapturederror", (ev) => console.error("WebGPU error:", (ev as GPUUncapturedErrorEvent).error.message));
494+ const g = makeGrid(params);
495+ $("info").innerHTML = `<div class="row"><span>Grid (default)</span><span>${g.nx} × ${g.ny} nodes</span></div><p class="sub">Press Start to run.</p>`;
496+ requestAnimationFrame(frame);
497+ applyUrlParams();
498+}
499+
500+/** URL parameters: any PhysParams key in SI units, plus autostart=1. Exposes the simulation on window for tests. */
501+function applyUrlParams(): void {
502+ const q = new URLSearchParams(location.search);
503+ for (const [k, v] of q) {
504+ if (k in params && k !== "radiation") {
505+ const num = parseFloat(v);
506+ if (isFinite(num)) {
507+ (params as unknown as Record<string, number>)[k] = num;
508+ const inp = inputs.get(k as keyof PhysParams);
509+ const spec = groups.flatMap((g) => g.fields).find((f) => f.key === k);
510+ if (inp && spec) (inp as HTMLInputElement).value = fmt(num / spec.scale);
511+ }
512+ }
513+ if (k === "radiation" && (v === "rayleigh" || v === "dipole")) params.radiation = v;
514+ if (k === "duration") {
515+ duration = params.duration;
516+ durInput.value = String(duration);
517+ }
518+ }
519+ updateDampingNote();
520+ (window as unknown as { steelsim: unknown }).steelsim = {
521+ status: () => (sim ? { ...sim.status(), dt: sim.dt, stepsPerSample: sim.stepsPerSample, nx: sim.grid.nx, ny: sim.grid.ny, stepsPerSec } : null),
522+ audio: () => (sim ? Array.from(sim.processedAudio()) : null),
523+ frame: (k: number) => (sim ? Array.from(sim.frames[k] ?? []) : null),
524+ frames: () => sim?.frames.length ?? 0,
525+ state: async () => { if (!sim) return null; const st = await sim.readState(); return { pos: Array.from(st.pos), vel: Array.from(st.vel) }; },
526+ pause: () => { paused = true; },
527+ canvasInfo: () => { const c = $<HTMLCanvasElement>("view"); return { cw: c.clientWidth, ch: c.clientHeight, w: c.width, h: c.height, mode: viewMode.value }; },
528+ setView: (m: ViewMode) => { viewMode.value = m; },
529+ scrub: (k: number) => { live = false; scrubFrame = k; scrub.value = String(k); updateLiveButton(); },
530+ };
531+ if (q.get("autostart") === "1") void start();
532+}
533+
534+function showNoGpu(msg: string): void {
535+ const el = $("nogpu");
536+ el.hidden = false;
537+ el.textContent = msg;
538+ btnStart.disabled = true;
539+}
540+
541+void init();
src/params.tsadded+114−0View file
@@ -0,0 +1,114 @@
1+// Physical and numerical parameters of the simulation, plus derived quantities.
2+
3+export interface PhysParams {
4+ Lx: number; // plate width along x [m]
5+ Ly: number; // plate height along y [m]
6+ h: number; // thickness [m]
7+ E: number; // Young's modulus [Pa]
8+ nu: number; // Poisson ratio
9+ rho: number; // density [kg/m^3]
10+ sigma0: number; // frequency-independent loss [1/s]
11+ sigma1: number; // frequency-dependent loss [m^2/s]
12+ holeInset: number; // distance of the thread hole from the two edges of the top-left corner [m]
13+ pushAmount: number; // corner displacement [m]
14+ fingerRadius: number; // radius of the pushed patch [m]
15+ rampTime: number; // duration of the push ramp [s]
16+ holdMax: number; // maximum settle time before release [s]
17+ settleSigma0: number; // artificial damping during push and settle [1/s]
18+ settleKE: number; // kinetic energy threshold that counts as settled [J]
19+ spacing: number; // requested grid spacing [m]
20+ fs: number; // audio sample rate [Hz]
21+ fps: number; // movie frame rate (also the batch rate) [Hz]
22+ duration: number; // recorded duration after release [s]
23+ listenerDist: number; // [m]
24+ listenerAzimuth: number; // [deg], rotation about the plate's vertical axis (0 = on the normal)
25+ listenerElevation: number; // [deg]
26+ radiation: "rayleigh" | "dipole";
27+}
28+
29+export const defaultParams: PhysParams = {
30+ Lx: 18 * 0.0254,
31+ Ly: 12 * 0.0254,
32+ h: 0.79e-3,
33+ E: 193e9,
34+ nu: 0.29,
35+ rho: 7900,
36+ sigma0: 1.15,
37+ sigma1: 1.08e-4,
38+ holeInset: 0.012,
39+ pushAmount: 0.01,
40+ fingerRadius: 0.01,
41+ rampTime: 0.1,
42+ holdMax: 0.5,
43+ settleSigma0: 300,
44+ settleKE: 1e-8,
45+ spacing: 0.004,
46+ fs: 48000,
47+ fps: 120,
48+ duration: 4,
49+ listenerDist: 1.0,
50+ listenerAzimuth: 0,
51+ listenerElevation: 0,
52+ radiation: "rayleigh",
53+};
54+
55+export const SPEED_OF_SOUND = 343; // [m/s]
56+export const RHO_AIR = 1.2; // [kg/m^3]
57+export const HIST_LEN = 256; // ring buffer length for the retarded-time audio sum [samples]
58+export const RB_INTERVAL = 32; // steps between rigid-body projections while ringing
59+export const RB_INTERVAL_HELD = 4; // while the finger constraint is active
60+
61+export interface Grid {
62+ nx: number;
63+ ny: number;
64+ hx: number;
65+ hy: number;
66+ n: number;
67+}
68+
69+export interface Material {
70+ D: number; // bending stiffness [N m]
71+ C: number; // membrane stiffness E h / (1 - nu^2) [N/m]
72+ G: number; // membrane shear stiffness E h / (2 (1 + nu)) [N/m]
73+ nu: number;
74+ rhoh: number; // mass per unit area [kg/m^2]
75+}
76+
77+export function makeGrid(p: PhysParams): Grid {
78+ const nx = Math.max(6, Math.round(p.Lx / p.spacing) + 1);
79+ const ny = Math.max(6, Math.round(p.Ly / p.spacing) + 1);
80+ return { nx, ny, hx: p.Lx / (nx - 1), hy: p.Ly / (ny - 1), n: nx * ny };
81+}
82+
83+export function makeMaterial(p: PhysParams): Material {
84+ return {
85+ D: (p.E * p.h ** 3) / (12 * (1 - p.nu * p.nu)),
86+ C: (p.E * p.h) / (1 - p.nu * p.nu),
87+ G: (p.E * p.h) / (2 * (1 + p.nu)),
88+ nu: p.nu,
89+ rhoh: p.rho * p.h,
90+ };
91+}
92+
93+/** Node coordinates relative to the plate center. */
94+export function nodeX(g: Grid, i: number): number {
95+ return (i - (g.nx - 1) / 2) * g.hx;
96+}
97+export function nodeY(g: Grid, j: number): number {
98+ return (j - (g.ny - 1) / 2) * g.hy;
99+}
100+
101+/** Trapezoidal weight of a node: 1 in the interior, 1/2 on an edge, 1/4 at a corner. */
102+export function trapWeight(g: Grid, i: number, j: number): number {
103+ const tx = i === 0 || i === g.nx - 1 ? 0.5 : 1;
104+ const ty = j === 0 || j === g.ny - 1 ? 0.5 : 1;
105+ return tx * ty;
106+}
107+
108+/** Decay time T60 [s] of a mode with angular frequency omega under the two-parameter loss model. */
109+export function t60(m: Material, p: PhysParams, omega: number): number {
110+ const kappa = Math.sqrt(m.D / m.rhoh);
111+ const k2 = omega / kappa;
112+ const sigma = p.sigma0 + p.sigma1 * k2;
113+ return (3 * Math.LN10) / sigma; // amplitude factor 1000 = 60 dB
114+}
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+168−0View file
@@ -0,0 +1,168 @@
1+// WebGPU rendering of the displacement field: a 2D color map or a lit 3D height field.
2+
3+import { type Grid, type PhysParams } 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+ camera = { yaw: 0.6, pitch: 0.5, dist: 1.0 };
16+
17+ private device: GPUDevice;
18+ private canvas: HTMLCanvasElement;
19+ private ctx: GPUCanvasContext;
20+ private format: GPUTextureFormat;
21+ private uniform: GPUBuffer;
22+ private pipe2d: GPURenderPipeline;
23+ private pipe3d: GPURenderPipeline;
24+ private bgLayout: GPUBindGroupLayout;
25+ private index: GPUBuffer;
26+ private indexCount: number;
27+ private depth: GPUTexture | null = null;
28+ private grid: Grid;
29+ private params: PhysParams;
30+ private bindGroups = new Map<GPUBuffer, GPUBindGroup>();
31+
32+ constructor(device: GPUDevice, canvas: HTMLCanvasElement, grid: Grid, params: PhysParams) {
33+ this.device = device;
34+ this.canvas = canvas;
35+ this.grid = grid;
36+ this.params = params;
37+ this.ctx = canvas.getContext("webgpu") as GPUCanvasContext;
38+ this.format = navigator.gpu.getPreferredCanvasFormat();
39+ this.ctx.configure({ device, format: this.format, alphaMode: "opaque" });
40+ this.uniform = device.createBuffer({ size: 128, usage: GPUBufferUsage.UNIFORM | GPUBufferUsage.COPY_DST });
41+ const module = device.createShaderModule({ code: shaderCode });
42+ this.bgLayout = device.createBindGroupLayout({
43+ entries: [
44+ { binding: 0, visibility: GPUShaderStage.VERTEX | GPUShaderStage.FRAGMENT, buffer: { type: "uniform" } },
45+ { binding: 1, visibility: GPUShaderStage.VERTEX | GPUShaderStage.FRAGMENT, buffer: { type: "read-only-storage" } },
46+ ],
47+ });
48+ const layout = device.createPipelineLayout({ bindGroupLayouts: [this.bgLayout] });
49+ this.pipe2d = device.createRenderPipeline({
50+ layout,
51+ vertex: { module, entryPoint: "vs2d" },
52+ fragment: { module, entryPoint: "fs2d", targets: [{ format: this.format }] },
53+ primitive: { topology: "triangle-list" },
54+ });
55+ this.pipe3d = device.createRenderPipeline({
56+ layout,
57+ vertex: { module, entryPoint: "vs3d" },
58+ fragment: { module, entryPoint: "fs3d", targets: [{ format: this.format }] },
59+ primitive: { topology: "triangle-list", cullMode: "none" },
60+ depthStencil: { format: "depth24plus", depthWriteEnabled: true, depthCompare: "less" },
61+ });
62+ // Triangle indices of the height field.
63+ const { nx, ny } = grid;
64+ const idx = new Uint32Array((nx - 1) * (ny - 1) * 6);
65+ let k = 0;
66+ for (let j = 0; j < ny - 1; j++) {
67+ for (let i = 0; i < nx - 1; i++) {
68+ const a = j * nx + i, b = a + 1, c = a + nx, d = c + 1;
69+ idx[k++] = a; idx[k++] = b; idx[k++] = d;
70+ idx[k++] = a; idx[k++] = d; idx[k++] = c;
71+ }
72+ }
73+ this.indexCount = idx.length;
74+ this.index = device.createBuffer({ size: idx.byteLength, usage: GPUBufferUsage.INDEX | GPUBufferUsage.COPY_DST });
75+ device.queue.writeBuffer(this.index, 0, idx);
76+ this.installControls();
77+ }
78+
79+ private installControls(): void {
80+ let dragging = false, lastX = 0, lastY = 0;
81+ this.canvas.addEventListener("pointerdown", (e) => {
82+ dragging = true; lastX = e.clientX; lastY = e.clientY;
83+ this.canvas.setPointerCapture(e.pointerId);
84+ });
85+ this.canvas.addEventListener("pointerup", () => (dragging = false));
86+ this.canvas.addEventListener("pointermove", (e) => {
87+ if (!dragging || this.mode !== "3d") return;
88+ this.camera.yaw += (e.clientX - lastX) * 0.01;
89+ this.camera.pitch = Math.max(-1.5, Math.min(1.5, this.camera.pitch + (e.clientY - lastY) * 0.01));
90+ lastX = e.clientX; lastY = e.clientY;
91+ });
92+ this.canvas.addEventListener("wheel", (e) => {
93+ if (this.mode !== "3d") return;
94+ e.preventDefault();
95+ this.camera.dist *= Math.exp(e.deltaY * 0.001);
96+ this.camera.dist = Math.max(0.2, Math.min(5, this.camera.dist));
97+ }, { passive: false });
98+ }
99+
100+ private resize(): void {
101+ const dpr = window.devicePixelRatio || 1;
102+ const w = Math.max(1, Math.floor(this.canvas.clientWidth * dpr));
103+ const h = Math.max(1, Math.floor(this.canvas.clientHeight * dpr));
104+ if (this.canvas.width !== w || this.canvas.height !== h || !this.depth) {
105+ this.canvas.width = w;
106+ this.canvas.height = h;
107+ this.depth?.destroy();
108+ this.depth = this.device.createTexture({ size: [w, h], format: "depth24plus", usage: GPUTextureUsage.RENDER_ATTACHMENT });
109+ }
110+ }
111+
112+ /** Renders the w field held in `wBuffer` (f32 per node). `wmax` drives the automatic color range. */
113+ render(encoder: GPUCommandEncoder, wBuffer: GPUBuffer, wmax: number): void {
114+ this.resize();
115+ const aspect = this.canvas.width / this.canvas.height;
116+ const range = this.colorRange > 0 ? this.colorRange : Math.max(wmax, 1e-6);
117+ const u = new Float32Array(32);
118+ const ui = new Uint32Array(u.buffer);
119+ const eye = [
120+ this.camera.dist * Math.cos(this.camera.pitch) * Math.sin(this.camera.yaw),
121+ this.camera.dist * Math.sin(this.camera.pitch),
122+ this.camera.dist * Math.cos(this.camera.pitch) * Math.cos(this.camera.yaw),
123+ ];
124+ const mvp = multiply(perspective(0.7, aspect, 0.02, 20), lookAt(eye, [0, 0, 0], [0, 1, 0]));
125+ u.set(mvp, 0);
126+ ui[16] = this.grid.nx; ui[17] = this.grid.ny;
127+ u[18] = this.grid.hx; u[19] = this.grid.hy;
128+ u[20] = this.exaggeration; u[21] = range; u[22] = aspect;
129+ ui[23] = this.mode === "3d" ? 1 : 0;
130+ const p = this.params;
131+ u[24] = -p.Lx / 2 + p.holeInset; u[25] = p.Ly / 2 - p.holeInset;
132+ u[26] = p.Lx / 2; u[27] = -p.Ly / 2;
133+ u[28] = p.fingerRadius;
134+ u[29] = 0.4; u[30] = 0.5; u[31] = 1.0;
135+ this.device.queue.writeBuffer(this.uniform, 0, u);
136+
137+ let bg = this.bindGroups.get(wBuffer);
138+ if (!bg) {
139+ bg = this.device.createBindGroup({
140+ layout: this.bgLayout,
141+ entries: [
142+ { binding: 0, resource: { buffer: this.uniform } },
143+ { binding: 1, resource: { buffer: wBuffer } },
144+ ],
145+ });
146+ this.bindGroups.set(wBuffer, bg);
147+ }
148+ const view = this.ctx.getCurrentTexture().createView();
149+ const clear = { r: 0.11, g: 0.12, b: 0.14, a: 1 };
150+ if (this.mode === "2d") {
151+ const pass = encoder.beginRenderPass({ colorAttachments: [{ view, clearValue: clear, loadOp: "clear", storeOp: "store" }] });
152+ pass.setPipeline(this.pipe2d);
153+ pass.setBindGroup(0, bg);
154+ pass.draw(3);
155+ pass.end();
156+ } else {
157+ const pass = encoder.beginRenderPass({
158+ colorAttachments: [{ view, clearValue: clear, loadOp: "clear", storeOp: "store" }],
159+ depthStencilAttachment: { view: this.depth!.createView(), depthClearValue: 1, depthLoadOp: "clear", depthStoreOp: "store" },
160+ });
161+ pass.setPipeline(this.pipe3d);
162+ pass.setBindGroup(0, bg);
163+ pass.setIndexBuffer(this.index, "uint32");
164+ pass.drawIndexed(this.indexCount);
165+ pass.end();
166+ }
167+ }
168+}
src/render/shaders.wgsladded+108−0View file
@@ -0,0 +1,108 @@
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+};
20+
21+@group(0) @binding(0) var<uniform> R: RU;
22+@group(0) @binding(1) var<storage, read> w: array<f32>;
23+
24+fn colormap(t: f32) -> vec3f {
25+ let s = clamp(t, -1.0, 1.0);
26+ let mid = vec3f(0.93, 0.93, 0.93);
27+ let blue = vec3f(0.10, 0.35, 0.85);
28+ let red = vec3f(0.85, 0.20, 0.12);
29+ if (s < 0.0) { return mix(mid, blue, -s); }
30+ return mix(mid, red, s);
31+}
32+
33+fn wAt(i: i32, j: i32) -> f32 {
34+ let ii = clamp(i, 0, i32(R.nx) - 1);
35+ let jj = clamp(j, 0, i32(R.ny) - 1);
36+ return w[u32(jj) * R.nx + u32(ii)];
37+}
38+
39+// ---------- 2D view: full-screen triangle, bilinear lookup in the fragment shader ----------
40+
41+struct V2 { @builtin(position) pos: vec4f, @location(0) uv: vec2f };
42+
43+@vertex
44+fn vs2d(@builtin(vertex_index) vi: u32) -> V2 {
45+ var p = array<vec2f, 3>(vec2f(-1.0, -1.0), vec2f(3.0, -1.0), vec2f(-1.0, 3.0));
46+ var o: V2;
47+ o.pos = vec4f(p[vi], 0.0, 1.0);
48+ o.uv = p[vi];
49+ return o;
50+}
51+
52+@fragment
53+fn fs2d(in: V2) -> @location(0) vec4f {
54+ let Lx = f32(R.nx - 1u) * R.hx;
55+ let Ly = f32(R.ny - 1u) * R.hy;
56+ // Fit the plate into the canvas with a margin, preserving its aspect ratio.
57+ let margin = 0.92;
58+ var sx = margin; var sy = margin;
59+ let plateAspect = Lx / Ly;
60+ if (plateAspect > R.aspect) { sy = margin * R.aspect / plateAspect; } else { sx = margin * plateAspect / R.aspect; }
61+ let px = in.uv.x / sx * 0.5 * Lx;
62+ let py = in.uv.y / sy * 0.5 * Ly;
63+ let bg = vec4f(0.11, 0.12, 0.14, 1.0);
64+ if (abs(px) > 0.5 * Lx || abs(py) > 0.5 * Ly) { return bg; }
65+ let fi = (px + 0.5 * Lx) / R.hx;
66+ let fj = (py + 0.5 * Ly) / R.hy;
67+ let i0 = i32(floor(fi)); let j0 = i32(floor(fj));
68+ let a = fi - f32(i0); let b = fj - f32(j0);
69+ let val = mix(mix(wAt(i0, j0), wAt(i0 + 1, j0), a), mix(wAt(i0, j0 + 1), wAt(i0 + 1, j0 + 1), a), b);
70+ var col = colormap(val / R.colorRange);
71+ // Markers: hole (ring) and pushed patch (dashed-looking ring).
72+ let dh = distance(vec2f(px, py), vec2f(R.holeX, R.holeY));
73+ if (abs(dh - R.markerR) < 0.35 * R.markerR) { col = vec3f(0.1, 0.1, 0.1); }
74+ if (dh < 0.4 * R.markerR) { col = vec3f(0.11, 0.12, 0.14); }
75+ let dp = distance(vec2f(px, py), vec2f(R.pushX, R.pushY));
76+ if (abs(dp - R.markerR * 1.6) < 0.25 * R.markerR) { col = vec3f(0.95, 0.75, 0.1); }
77+ return vec4f(col, 1.0);
78+}
79+
80+// ---------- 3D view: height field ----------
81+
82+struct V3 { @builtin(position) pos: vec4f, @location(0) val: f32, @location(1) normal: vec3f };
83+
84+@vertex
85+fn vs3d(@builtin(vertex_index) vi: u32) -> V3 {
86+ let i = i32(vi % R.nx); let j = i32(vi / R.nx);
87+ let x = (f32(i) - 0.5 * f32(R.nx - 1u)) * R.hx;
88+ let y = (f32(j) - 0.5 * f32(R.ny - 1u)) * R.hy;
89+ let val = wAt(i, j);
90+ let z = val * R.wScale;
91+ let dzdx = (wAt(i + 1, j) - wAt(i - 1, j)) * R.wScale / (2.0 * R.hx);
92+ let dzdy = (wAt(i, j + 1) - wAt(i, j - 1)) * R.wScale / (2.0 * R.hy);
93+ var o: V3;
94+ o.pos = R.mvp * vec4f(x, y, z, 1.0);
95+ o.val = val;
96+ o.normal = normalize(vec3f(-dzdx, -dzdy, 1.0));
97+ return o;
98+}
99+
100+@fragment
101+fn fs3d(in: V3) -> @location(0) vec4f {
102+ let l = normalize(vec3f(R.lightX, R.lightY, R.lightZ));
103+ let n = normalize(in.normal);
104+ let diff = abs(dot(n, l));
105+ let base = colormap(in.val / R.colorRange);
106+ let col = base * (0.35 + 0.65 * diff) + vec3f(0.15) * pow(diff, 24.0);
107+ return vec4f(col, 1.0);
108+}
src/sim/cpuRef.tsadded+199−0View file
@@ -0,0 +1,199 @@
1+// CPU reference implementation of the spatial discretization (float64).
2+// It mirrors step.wgsl exactly and is used for the stable time-step estimate
3+// and for the tests. It is not used in the time loop of the browser simulation.
4+
5+import { type Grid, type Material, trapWeight } from "../params.ts";
6+
7+export interface Fields {
8+ u: Float64Array;
9+ v: Float64Array;
10+ w: Float64Array;
11+}
12+
13+const G3 = 1 / Math.sqrt(3);
14+
15+/**
16+ * Accelerations (au, av, aw) of every node for displacement fields (u, v, w),
17+ * without damping. Bending forces come from the discrete Kirchhoff energy
18+ * (natural free-edge conditions); membrane forces come from bilinear (Q4)
19+ * elements with 2x2 Gauss integration and the von Karman strains.
20+ */
21+export function accelerations(g: Grid, m: Material, f: Fields, out: Fields): void {
22+ const { nx, ny, hx, hy } = g;
23+ const { u, v, w } = f;
24+ const hx2 = hx * hx;
25+ const hy2 = hy * hy;
26+ const nu = m.nu;
27+ const id = (i: number, j: number) => j * nx + i;
28+
29+ // Bending moment-like quantities at nodes (zero where the second difference is undefined).
30+ const Ma = (i: number, j: number): number => {
31+ if (i < 1 || i > nx - 2 || j < 0 || j > ny - 1) return 0;
32+ const a = (w[id(i + 1, j)] - 2 * w[id(i, j)] + w[id(i - 1, j)]) / hx2;
33+ let b = 0;
34+ if (j >= 1 && j <= ny - 2) b = (w[id(i, j + 1)] - 2 * w[id(i, j)] + w[id(i, j - 1)]) / hy2;
35+ return trapWeight(g, i, j) * (a + nu * b);
36+ };
37+ const Mb = (i: number, j: number): number => {
38+ if (j < 1 || j > ny - 2 || i < 0 || i > nx - 1) return 0;
39+ const b = (w[id(i, j + 1)] - 2 * w[id(i, j)] + w[id(i, j - 1)]) / hy2;
40+ let a = 0;
41+ if (i >= 1 && i <= nx - 2) a = (w[id(i + 1, j)] - 2 * w[id(i, j)] + w[id(i - 1, j)]) / hx2;
42+ return trapWeight(g, i, j) * (b + nu * a);
43+ };
44+ // Mixed second difference at a cell center (zero for cells outside the plate).
45+ const Cc = (ci: number, cj: number): number => {
46+ if (ci < 0 || ci > nx - 2 || cj < 0 || cj > ny - 2) return 0;
47+ return (w[id(ci + 1, cj + 1)] - w[id(ci + 1, cj)] - w[id(ci, cj + 1)] + w[id(ci, cj)]) / (hx * hy);
48+ };
49+
50+ const q = (hx * hy) / 4; // area per Gauss point
51+ const U = new Float64Array(4);
52+ const V = new Float64Array(4);
53+ const Wn = new Float64Array(4);
54+ const dNdx = new Float64Array(4);
55+ const dNdy = new Float64Array(4);
56+
57+ // Membrane force on local node `a` of cell (ci, cj); returns [Fu, Fv, Fw].
58+ const cellForce = (ci: number, cj: number, a: number, F: Float64Array): void => {
59+ if (ci < 0 || ci > nx - 2 || cj < 0 || cj > ny - 2) return;
60+ const n0 = id(ci, cj), n1 = id(ci + 1, cj), n2 = id(ci + 1, cj + 1), n3 = id(ci, cj + 1);
61+ U[0] = u[n0]; U[1] = u[n1]; U[2] = u[n2]; U[3] = u[n3];
62+ V[0] = v[n0]; V[1] = v[n1]; V[2] = v[n2]; V[3] = v[n3];
63+ Wn[0] = w[n0]; Wn[1] = w[n1]; Wn[2] = w[n2]; Wn[3] = w[n3];
64+ for (let gp = 0; gp < 4; gp++) {
65+ const xi = gp & 1 ? G3 : -G3;
66+ const eta = gp & 2 ? G3 : -G3;
67+ dNdx[0] = (-(1 - eta) / 4) * (2 / hx);
68+ dNdx[1] = ((1 - eta) / 4) * (2 / hx);
69+ dNdx[2] = ((1 + eta) / 4) * (2 / hx);
70+ dNdx[3] = (-(1 + eta) / 4) * (2 / hx);
71+ dNdy[0] = (-(1 - xi) / 4) * (2 / hy);
72+ dNdy[1] = (-(1 + xi) / 4) * (2 / hy);
73+ dNdy[2] = ((1 + xi) / 4) * (2 / hy);
74+ dNdy[3] = ((1 - xi) / 4) * (2 / hy);
75+ let ux = 0, uy = 0, vx = 0, vy = 0, wx = 0, wy = 0;
76+ for (let b = 0; b < 4; b++) {
77+ ux += dNdx[b] * U[b]; uy += dNdy[b] * U[b];
78+ vx += dNdx[b] * V[b]; vy += dNdy[b] * V[b];
79+ wx += dNdx[b] * Wn[b]; wy += dNdy[b] * Wn[b];
80+ }
81+ const exx = ux + 0.5 * wx * wx;
82+ const eyy = vy + 0.5 * wy * wy;
83+ const gxy = uy + vx + wx * wy;
84+ const Nxx = m.C * (exx + nu * eyy);
85+ const Nyy = m.C * (eyy + nu * exx);
86+ const Nxy = m.G * gxy;
87+ F[0] -= q * (Nxx * dNdx[a] + Nxy * dNdy[a]);
88+ F[1] -= q * (Nyy * dNdy[a] + Nxy * dNdx[a]);
89+ F[2] -= q * ((Nxx * wx + Nxy * wy) * dNdx[a] + (Nxy * wx + Nyy * wy) * dNdy[a]);
90+ }
91+ };
92+
93+ const F = new Float64Array(3);
94+ for (let j = 0; j < ny; j++) {
95+ for (let i = 0; i < nx; i++) {
96+ const k = id(i, j);
97+ const bend =
98+ (Ma(i - 1, j) - 2 * Ma(i, j) + Ma(i + 1, j)) / hx2 +
99+ (Mb(i, j - 1) - 2 * Mb(i, j) + Mb(i, j + 1)) / hy2 +
100+ (2 * (1 - nu) * (Cc(i, j) - Cc(i - 1, j) - Cc(i, j - 1) + Cc(i - 1, j - 1))) / (hx * hy);
101+ F[0] = 0; F[1] = 0; F[2] = -m.D * hx * hy * bend;
102+ cellForce(i, j, 0, F);
103+ cellForce(i - 1, j, 1, F);
104+ cellForce(i - 1, j - 1, 2, F);
105+ cellForce(i, j - 1, 3, F);
106+ const mass = m.rhoh * hx * hy * trapWeight(g, i, j);
107+ out.u[k] = F[0] / mass;
108+ out.v[k] = F[1] / mass;
109+ out.w[k] = F[2] / mass;
110+ }
111+ }
112+}
113+
114+/**
115+ * Largest angular frequency of the linearized discrete system, by power
116+ * iteration on a small grid with the same spacing (the extreme mode is a
117+ * grid-scale mode, so the plate size does not matter).
118+ */
119+export function maxAngularFrequency(g: Grid, m: Material, nSmall = 12, iters = 300): number {
120+ const gs: Grid = { nx: nSmall, ny: nSmall, hx: g.hx, hy: g.hy, n: nSmall * nSmall };
121+ const n = gs.n;
122+ const x: Fields = { u: new Float64Array(n), v: new Float64Array(n), w: new Float64Array(n) };
123+ const y: Fields = { u: new Float64Array(n), v: new Float64Array(n), w: new Float64Array(n) };
124+ let seed = 12345;
125+ const rnd = () => {
126+ seed = (seed * 1664525 + 1013904223) >>> 0;
127+ return seed / 4294967296 - 0.5;
128+ };
129+ for (let k = 0; k < n; k++) {
130+ x.u[k] = rnd(); x.v[k] = rnd(); x.w[k] = rnd();
131+ }
132+ const eps = 1e-9; // keeps the nonlinear terms negligible
133+ let lambda = 0;
134+ for (let it = 0; it < iters; it++) {
135+ let norm = 0;
136+ for (let k = 0; k < n; k++) norm += x.u[k] ** 2 + x.v[k] ** 2 + x.w[k] ** 2;
137+ norm = Math.sqrt(norm);
138+ for (let k = 0; k < n; k++) {
139+ x.u[k] *= eps / norm; x.v[k] *= eps / norm; x.w[k] *= eps / norm;
140+ }
141+ accelerations(gs, m, x, y);
142+ let dot = 0;
143+ 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];
144+ lambda = -dot / (eps * eps); // Rayleigh quotient of -M^{-1}K
145+ for (let k = 0; k < n; k++) {
146+ x.u[k] = -y.u[k]; x.v[k] = -y.v[k]; x.w[k] = -y.w[k];
147+ }
148+ }
149+ return Math.sqrt(Math.max(lambda, 0));
150+}
151+
152+/** Discrete potential energy (bending + membrane) matching `accelerations`. */
153+export function potentialEnergy(g: Grid, m: Material, f: Fields): number {
154+ const { nx, ny, hx, hy } = g;
155+ const { u, v, w } = f;
156+ const nu = m.nu;
157+ const id = (i: number, j: number) => j * nx + i;
158+ let Ub = 0;
159+ for (let j = 0; j < ny; j++) {
160+ for (let i = 0; i < nx; i++) {
161+ let a = 0, b = 0;
162+ if (i >= 1 && i <= nx - 2) a = (w[id(i + 1, j)] - 2 * w[id(i, j)] + w[id(i - 1, j)]) / (hx * hx);
163+ if (j >= 1 && j <= ny - 2) b = (w[id(i, j + 1)] - 2 * w[id(i, j)] + w[id(i, j - 1)]) / (hy * hy);
164+ Ub += trapWeight(g, i, j) * (a * a + b * b + 2 * nu * a * b);
165+ }
166+ }
167+ for (let cj = 0; cj < ny - 1; cj++) {
168+ for (let ci = 0; ci < nx - 1; ci++) {
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+ for (let cj = 0; cj < ny - 1; cj++) {
178+ for (let ci = 0; ci < nx - 1; ci++) {
179+ const nn = [id(ci, cj), id(ci + 1, cj), id(ci + 1, cj + 1), id(ci, cj + 1)];
180+ for (let gp = 0; gp < 4; gp++) {
181+ const xi = gp & 1 ? G3 : -G3;
182+ const eta = gp & 2 ? G3 : -G3;
183+ const dNdx = [-(1 - eta), 1 - eta, 1 + eta, -(1 + eta)].map((s) => (s / 4) * (2 / hx));
184+ const dNdy = [-(1 - xi), -(1 + xi), 1 + xi, 1 - xi].map((s) => (s / 4) * (2 / hy));
185+ let ux = 0, uy = 0, vx = 0, vy = 0, wx = 0, wy = 0;
186+ for (let b = 0; b < 4; b++) {
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]]; wy += dNdy[b] * w[nn[b]];
190+ }
191+ const exx = ux + 0.5 * wx * wx;
192+ const eyy = vy + 0.5 * wy * wy;
193+ const gxy = uy + vx + wx * wy;
194+ Um += q * 0.5 * (m.C * (exx * exx + eyy * eyy + 2 * nu * exx * eyy) + m.G * gxy * gxy);
195+ }
196+ }
197+ }
198+ return Ub + Um;
199+}
src/sim/shaders/audio.wgsladded+36−0View file
@@ -0,0 +1,36 @@
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+
5+@group(0) @binding(0) var<uniform> P: Params;
6+@group(0) @binding(1) var<uniform> SP: StepParams;
7+@group(0) @binding(2) var<storage, read> hist: array<f32>;
8+@group(0) @binding(3) var<storage, read> lw: array<vec2f>;
9+@group(0) @binding(4) var<storage, read_write> audioOut: array<f32>;
10+
11+var<workgroup> sh: array<f32, 256>;
12+
13+@compute @workgroup_size(256)
14+fn main(@builtin(local_invocation_id) lid: vec3u) {
15+ let t = lid.x;
16+ let H = i32(P.histLen);
17+ let slot = SP.audioSlot;
18+ var acc = 0.0;
19+ for (var k = t; k < P.n; k += 256u) {
20+ let d = lw[k].y;
21+ let k0 = i32(floor(d));
22+ let f = d - f32(k0);
23+ let s0 = ((slot - k0) % H + H) % H;
24+ let s1 = ((slot - k0 - 1) % H + H) % H;
25+ let a0 = hist[u32(s0) * P.n + k];
26+ let a1 = hist[u32(s1) * P.n + k];
27+ acc += lw[k].x * mix(a0, a1, f);
28+ }
29+ sh[t] = acc;
30+ workgroupBarrier();
31+ for (var stride = 128u; stride > 0u; stride >>= 1u) {
32+ if (t < stride) { sh[t] += sh[t + stride]; }
33+ workgroupBarrier();
34+ }
35+ if (t == 0u) { audioOut[u32(SP.audioOut)] = sh[0]; }
36+}
src/sim/shaders/common.wgsladded+45−0View file
@@ -0,0 +1,45 @@
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+ pad0: f32,
24+};
25+
26+// Per-step parameters, addressed with a dynamic uniform offset.
27+struct StepParams {
28+ step: u32,
29+ flags: u32, // bit 0: push constraint active
30+ audioSlot: i32, // ring-buffer slot to record acceleration into (-1: none)
31+ audioOut: i32, // output sample index within the batch
32+ pushW: f32, // prescribed patch displacement [m]
33+ pushV: f32, // prescribed patch velocity [m/s]
34+ sigma0: f32,
35+ sigma1: f32,
36+};
37+
38+fn nodeX(i: i32) -> f32 { return (f32(i) - 0.5 * f32(P.nx - 1u)) * P.hx; }
39+fn nodeY(j: i32) -> f32 { return (f32(j) - 0.5 * f32(P.ny - 1u)) * P.hy; }
40+fn trapWeight(i: i32, j: i32) -> f32 {
41+ var t = 1.0;
42+ if (i == 0 || i == i32(P.nx) - 1) { t = t * 0.5; }
43+ if (j == 0 || j == i32(P.ny) - 1) { t = t * 0.5; }
44+ return t;
45+}
src/sim/shaders/extract.wgsladded+12−0View file
@@ -0,0 +1,12 @@
1+// Copies the out-of-plane 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+68−0View file
@@ -0,0 +1,68 @@
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 * trapWeight(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+ let x = nodeX(i); let y = nodeY(j);
61+ let n = P.n;
62+ pos[k] -= stats[0] - stats[3] * y;
63+ pos[n + k] -= stats[1] + stats[3] * x;
64+ pos[2u * n + k] -= stats[2] + stats[4] * y - stats[5] * x;
65+ vel[k] -= stats[6] - stats[9] * y;
66+ vel[n + k] -= stats[7] + stats[9] * x;
67+ vel[2u * n + k] -= stats[8] + stats[10] * y - stats[11] * x;
68+}
src/sim/shaders/step.wgsladded+175−0View file
@@ -0,0 +1,175 @@
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 u, v, w, dw into
7+// shared memory once, so the stencils below never touch global memory.
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> posIn: array<f32>;
12+@group(0) @binding(3) var<storage, read> velIn: array<f32>;
13+@group(0) @binding(4) var<storage, read_write> posOut: array<f32>;
14+@group(0) @binding(5) var<storage, read_write> velOut: array<f32>;
15+@group(0) @binding(6) var<storage, read_write> hist: array<f32>; // [slot * n + node] acceleration of w
16+
17+const TILE = 16;
18+const HALO = 2;
19+const TW = 20; // TILE + 2 * HALO
20+var<workgroup> tw: array<f32, 400>;
21+// Membrane force contributions of the 17x17 cells around the tile, 4 nodes x 3 components each.
22+const CT = 17;
23+var<workgroup> cellF: array<f32, 3468>;
24+
25+// Tile origin in global node coordinates (set once per invocation).
26+var<private> ox: i32;
27+var<private> oy: i32;
28+var<private> ihx2: f32;
29+var<private> ihy2: f32;
30+var<private> ihxy: f32;
31+
32+fn id(i: i32, j: i32) -> u32 { return u32(j) * P.nx + u32(i); }
33+fn lidx(i: i32, j: i32) -> u32 { return u32(j - oy) * u32(TW) + u32(i - ox); }
34+// Global indices; the tile covers [ox, ox + TW) x [oy, oy + TW).
35+fn W(i: i32, j: i32) -> f32 { return tw[lidx(i, j)]; }
36+fn U(i: i32, j: i32) -> f32 { return posIn[id(i, j)]; }
37+fn V(i: i32, j: i32) -> f32 { return posIn[P.n + id(i, j)]; }
38+fn DW(i: i32, j: i32) -> f32 { return velIn[2u * P.n + id(i, j)]; }
39+
40+// Weighted bending "moments" at nodes; zero where the second difference is undefined.
41+fn Ma(i: i32, j: i32) -> f32 {
42+ let nx = i32(P.nx); let ny = i32(P.ny);
43+ if (i < 1 || i > nx - 2 || j < 0 || j > ny - 1) { return 0.0; }
44+ let a = (W(i + 1, j) - 2.0 * W(i, j) + W(i - 1, j)) * ihx2;
45+ var b = 0.0;
46+ if (j >= 1 && j <= ny - 2) { b = (W(i, j + 1) - 2.0 * W(i, j) + W(i, j - 1)) * ihy2; }
47+ return trapWeight(i, j) * (a + P.nu * b);
48+}
49+fn Mb(i: i32, j: i32) -> f32 {
50+ let nx = i32(P.nx); let ny = i32(P.ny);
51+ if (j < 1 || j > ny - 2 || i < 0 || i > nx - 1) { return 0.0; }
52+ let b = (W(i, j + 1) - 2.0 * W(i, j) + W(i, j - 1)) * ihy2;
53+ var a = 0.0;
54+ if (i >= 1 && i <= nx - 2) { a = (W(i + 1, j) - 2.0 * W(i, j) + W(i - 1, j)) * ihx2; }
55+ return trapWeight(i, j) * (b + P.nu * a);
56+}
57+// Mixed second difference at a cell center; zero outside the plate.
58+fn Cc(ci: i32, cj: i32) -> f32 {
59+ if (ci < 0 || ci > i32(P.nx) - 2 || cj < 0 || cj > i32(P.ny) - 2) { return 0.0; }
60+ return (W(ci + 1, cj + 1) - W(ci + 1, cj) - W(ci, cj + 1) + W(ci, cj)) * ihxy;
61+}
62+
63+// Membrane forces of the bilinear cell (ci, cj) on its four nodes (2x2 Gauss),
64+// written to cellF at local cell index lc.
65+fn cellForces(ci: i32, cj: i32, lc: u32) {
66+ let base = lc * 12u;
67+ if (ci < 0 || ci > i32(P.nx) - 2 || cj < 0 || cj > i32(P.ny) - 2) {
68+ for (var q = 0u; q < 12u; q++) { cellF[base + q] = 0.0; }
69+ return;
70+ }
71+ let Un = vec4f(U(ci, cj), U(ci + 1, cj), U(ci + 1, cj + 1), U(ci, cj + 1));
72+ let Vn = vec4f(V(ci, cj), V(ci + 1, cj), V(ci + 1, cj + 1), V(ci, cj + 1));
73+ let Wn = vec4f(W(ci, cj), W(ci + 1, cj), W(ci + 1, cj + 1), W(ci, cj + 1));
74+ let q = 0.25 * P.hx * P.hy;
75+ let ihx = 1.0 / P.hx; let ihy = 1.0 / P.hy;
76+ let g3 = 0.57735026918962576;
77+ var Fu = vec4f(0.0); var Fv = vec4f(0.0); var Fw = vec4f(0.0);
78+ for (var gp = 0; gp < 4; gp++) {
79+ var xi = -g3; if ((gp & 1) != 0) { xi = g3; }
80+ var eta = -g3; if ((gp & 2) != 0) { eta = g3; }
81+ let dNdx = vec4f(-(1.0 - eta), 1.0 - eta, 1.0 + eta, -(1.0 + eta)) * (0.5 * ihx);
82+ let dNdy = vec4f(-(1.0 - xi), -(1.0 + xi), 1.0 + xi, 1.0 - xi) * (0.5 * ihy);
83+ let ux = dot(dNdx, Un); let uy = dot(dNdy, Un);
84+ let vx = dot(dNdx, Vn); let vy = dot(dNdy, Vn);
85+ let wx = dot(dNdx, Wn); let wy = dot(dNdy, Wn);
86+ let exx = ux + 0.5 * wx * wx;
87+ let eyy = vy + 0.5 * wy * wy;
88+ let gxy = uy + vx + wx * wy;
89+ let Nxx = P.C * (exx + P.nu * eyy);
90+ let Nyy = P.C * (eyy + P.nu * exx);
91+ let Nxy = P.G * gxy;
92+ Fu -= q * (Nxx * dNdx + Nxy * dNdy);
93+ Fv -= q * (Nyy * dNdy + Nxy * dNdx);
94+ Fw -= q * ((Nxx * wx + Nxy * wy) * dNdx + (Nxy * wx + Nyy * wy) * dNdy);
95+ }
96+ for (var a = 0u; a < 4u; a++) {
97+ cellF[base + a * 3u] = Fu[a];
98+ cellF[base + a * 3u + 1u] = Fv[a];
99+ cellF[base + a * 3u + 2u] = Fw[a];
100+ }
101+}
102+
103+// Stored membrane force on local node `a` of cell (ci, cj).
104+fn cellForce(ci: i32, cj: i32, a: u32) -> vec3f {
105+ let lc = u32(cj - oy - 1) * u32(CT) + u32(ci - ox - 1);
106+ let base = lc * 12u + a * 3u;
107+ return vec3f(cellF[base], cellF[base + 1u], cellF[base + 2u]);
108+}
109+
110+@compute @workgroup_size(16, 16)
111+fn main(@builtin(global_invocation_id) gid: vec3u, @builtin(local_invocation_id) lid: vec3u, @builtin(workgroup_id) wid: vec3u) {
112+ let n = P.n;
113+ ox = i32(wid.x) * TILE - HALO;
114+ oy = i32(wid.y) * TILE - HALO;
115+ ihx2 = 1.0 / (P.hx * P.hx); ihy2 = 1.0 / (P.hy * P.hy); ihxy = 1.0 / (P.hx * P.hy);
116+ // Cooperative load of the w tile (400 entries, 256 threads).
117+ let lt = lid.y * 16u + lid.x;
118+ for (var e = lt; e < 400u; e += 256u) {
119+ let gi = ox + i32(e % u32(TW)); let gj = oy + i32(e / u32(TW));
120+ var ww = 0.0;
121+ if (gi >= 0 && gi < i32(P.nx) && gj >= 0 && gj < i32(P.ny)) { ww = posIn[2u * n + id(gi, gj)]; }
122+ tw[e] = ww;
123+ }
124+ workgroupBarrier();
125+ // Membrane forces of the 17x17 cells adjacent to the tile's nodes.
126+ for (var e = lt; e < 289u; e += 256u) {
127+ let lci = i32(e % u32(CT)); let lcj = i32(e / u32(CT));
128+ cellForces(ox + 1 + lci, oy + 1 + lcj, e);
129+ }
130+ workgroupBarrier();
131+
132+ let i = i32(gid.x); let j = i32(gid.y);
133+ if (i >= i32(P.nx) || j >= i32(P.ny)) { return; }
134+ let k = id(i, j);
135+
136+
137+ let bend = (Ma(i - 1, j) - 2.0 * Ma(i, j) + Ma(i + 1, j)) * ihx2
138+ + (Mb(i, j - 1) - 2.0 * Mb(i, j) + Mb(i, j + 1)) * ihy2
139+ + 2.0 * (1.0 - P.nu) * (Cc(i, j) - Cc(i - 1, j) - Cc(i, j - 1) + Cc(i - 1, j - 1)) * ihxy;
140+ var F = vec3f(0.0, 0.0, -P.D * P.hx * P.hy * bend);
141+ F += cellForce(i, j, 0u);
142+ F += cellForce(i - 1, j, 1u);
143+ F += cellForce(i - 1, j - 1, 2u);
144+ F += cellForce(i, j - 1, 3u);
145+ let mass = P.rhoh * P.hx * P.hy * trapWeight(i, j);
146+ var acc = F * (1.0 / mass);
147+
148+ // Loss: frequency-independent on all components, frequency-dependent (graph Laplacian) on w.
149+ let vel = vec3f(velIn[k], velIn[n + k], velIn[2u * n + k]);
150+ var lap = 0.0;
151+ if (i > 0) { lap += (DW(i - 1, j) - vel.z) * ihx2; }
152+ if (i < i32(P.nx) - 1) { lap += (DW(i + 1, j) - vel.z) * ihx2; }
153+ if (j > 0) { lap += (DW(i, j - 1) - vel.z) * ihy2; }
154+ if (j < i32(P.ny) - 1) { lap += (DW(i, j + 1) - vel.z) * ihy2; }
155+ acc -= 2.0 * SP.sigma0 * vel;
156+ acc.z += 2.0 * SP.sigma1 * lap;
157+
158+ if (SP.audioSlot >= 0) { hist[u32(SP.audioSlot) * n + k] = acc.z; }
159+
160+ var vNew = vel + P.dt * acc;
161+ let pOld = vec3f(posIn[k], posIn[n + k], posIn[2u * n + k]);
162+ var pNew = pOld + P.dt * vNew;
163+ // Compensated update of w: the increment dt * v can be far below the
164+ // float32 resolution of w, so the rounding error is carried separately.
165+ let y = P.dt * vNew.z + posIn[3u * n + k];
166+ let t = pOld.z + y;
167+ var lo = (pOld.z - t) + y;
168+ pNew.z = t;
169+ if ((SP.flags & 1u) != 0u) {
170+ let dxp = nodeX(i) - P.pushX; let dyp = nodeY(j) - P.pushY;
171+ if (dxp * dxp + dyp * dyp <= P.pushR2) { pNew.z = SP.pushW; vNew.z = SP.pushV; lo = 0.0; }
172+ }
173+ posOut[k] = pNew.x; posOut[n + k] = pNew.y; posOut[2u * n + k] = pNew.z; posOut[3u * n + k] = lo;
174+ velOut[k] = vNew.x; velOut[n + k] = vNew.y; velOut[2u * n + k] = vNew.z;
175+}
src/sim/simulation.tsadded+538−0View file
@@ -0,0 +1,538 @@
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 batch
4+// encodes one compute dispatch per time step (with a per-step uniform addressed
5+// by dynamic offset), a rigid-body projection every RB_INTERVAL steps, one
6+// retarded-time audio sum per audio sample, and one movie-frame extraction at
7+// the end. Audio, frame, and statistics are read back after every batch.
8+
9+import {
10+ type PhysParams,
11+ type Grid,
12+ type Material,
13+ makeGrid,
14+ makeMaterial,
15+ nodeX,
16+ nodeY,
17+ trapWeight,
18+ HIST_LEN,
19+ RB_INTERVAL,
20+ RB_INTERVAL_HELD,
21+ RHO_AIR,
22+ SPEED_OF_SOUND,
23+} from "../params.ts";
24+import { maxAngularFrequency } from "./cpuRef.ts";
25+import commonWgsl from "./shaders/common.wgsl?raw";
26+import stepWgsl from "./shaders/step.wgsl?raw";
27+import rigidWgsl from "./shaders/rigid.wgsl?raw";
28+import audioWgsl from "./shaders/audio.wgsl?raw";
29+import extractWgsl from "./shaders/extract.wgsl?raw";
30+
31+export type Phase = "push" | "hold" | "recording" | "done";
32+
33+export interface SimStatus {
34+ phase: Phase;
35+ /** Simulation time relative to the release [s]. */
36+ t: number;
37+ totalSteps: number;
38+ ke: number;
39+ wmax: number;
40+ samples: number;
41+ frames: number;
42+ peakPressure: number;
43+ delayClamped: boolean;
44+}
45+
46+const STATS_LEN = 16;
47+
48+export class Simulation {
49+ readonly grid: Grid;
50+ readonly mat: Material;
51+ readonly dt: number;
52+ /** Time steps per audio sample. */
53+ readonly stepsPerSample: number;
54+ readonly samplesPerBatch: number;
55+ readonly stepsPerBatch: number;
56+ readonly omegaMax: number;
57+
58+ /** Raw listener pressure [Pa] at fs, from the release. */
59+ audio: Float32Array = new Float32Array(0);
60+ samples = 0;
61+ /** Movie frames of w [m]; frame k is at t = k / fps. */
62+ frames: Float32Array<ArrayBuffer>[] = [];
63+
64+ phase: Phase = "push";
65+ /** Steps since the start of the current phase clock (prep or recording). */
66+ private stepCount = 0;
67+ private prepTime = 0;
68+ totalSteps = 0;
69+ ke = 0;
70+ wmax = 0;
71+ peakPressure = 0;
72+ delayClamped = false;
73+ private hpState = { y: 0, x: 0 };
74+
75+ private device: GPUDevice;
76+ private params: PhysParams;
77+ private paramsBuf: GPUBuffer;
78+ private stepBuf: GPUBuffer;
79+ private stepStride: number;
80+ private stepData: ArrayBuffer;
81+ private pos: GPUBuffer[];
82+ private vel: GPUBuffer[];
83+ private hist: GPUBuffer;
84+ private lw: GPUBuffer;
85+ private audioOut: GPUBuffer;
86+ private frameBuf: GPUBuffer;
87+ private statsBuf: GPUBuffer;
88+ private staging: GPUBuffer[];
89+ private stagingIdx = 0;
90+ /** Index of the buffer pair holding the current state. */
91+ private cur = 0;
92+ /** The `displayW` buffer is a plain f32 array of w for the renderer. */
93+ readonly displayW: GPUBuffer;
94+
95+ private stepPipe: GPUComputePipeline;
96+ private sumsPipe: GPUComputePipeline;
97+ private applyPipe: GPUComputePipeline;
98+ private audioPipe: GPUComputePipeline;
99+ private extractPipe: GPUComputePipeline;
100+ private stepBG: GPUBindGroup[];
101+ private rbBG: GPUBindGroup[];
102+ private audioBG: GPUBindGroup;
103+ private extractBG: GPUBindGroup[];
104+ private displayBG: GPUBindGroup[];
105+ private stagingSize: number;
106+
107+ constructor(device: GPUDevice, params: PhysParams) {
108+ this.device = device;
109+ this.params = { ...params };
110+ const g = makeGrid(params);
111+ const m = makeMaterial(params);
112+ this.grid = g;
113+ this.mat = m;
114+
115+ // Stable time step, aligned so that an integer number of steps make one audio sample.
116+ this.omegaMax = maxAngularFrequency(g, m);
117+ const dtMax = (0.85 * 2) / this.omegaMax;
118+ this.stepsPerSample = Math.ceil(1 / (params.fs * dtMax));
119+ this.dt = 1 / (params.fs * this.stepsPerSample);
120+ this.samplesPerBatch = Math.max(1, Math.round(params.fs / params.fps));
121+ this.stepsPerBatch = this.samplesPerBatch * this.stepsPerSample;
122+
123+ const n = g.n;
124+ const f32 = (size: number, usage = GPUBufferUsage.STORAGE | GPUBufferUsage.COPY_SRC | GPUBufferUsage.COPY_DST) =>
125+ device.createBuffer({ size: Math.max(16, size * 4), usage });
126+
127+ this.pos = [f32(4 * n), f32(4 * n)]; // [u | v | w | wLo]
128+ this.vel = [f32(3 * n), f32(3 * n)]; // [du | dv | dw]
129+ this.hist = f32(HIST_LEN * n);
130+ this.lw = f32(2 * n);
131+ this.audioOut = f32(this.samplesPerBatch);
132+ this.frameBuf = f32(n);
133+ this.displayW = f32(n);
134+ this.statsBuf = f32(STATS_LEN);
135+ this.stagingSize = (this.samplesPerBatch + n + STATS_LEN) * 4;
136+ this.staging = [0, 1].map(() =>
137+ device.createBuffer({ size: this.stagingSize, usage: GPUBufferUsage.MAP_READ | GPUBufferUsage.COPY_DST }),
138+ );
139+
140+ // Constant parameters.
141+ let mass = 0, Ix = 0, Iy = 0;
142+ for (let j = 0; j < g.ny; j++) {
143+ for (let i = 0; i < g.nx; i++) {
144+ const mk = m.rhoh * g.hx * g.hy * trapWeight(g, i, j);
145+ mass += mk;
146+ Ix += mk * nodeY(g, j) ** 2;
147+ Iy += mk * nodeX(g, i) ** 2;
148+ }
149+ }
150+ const pf = new Float32Array(20);
151+ const pu = new Uint32Array(pf.buffer);
152+ pu[0] = g.nx; pu[1] = g.ny; pu[2] = n; pu[3] = HIST_LEN;
153+ pf[4] = g.hx; pf[5] = g.hy; pf[6] = this.dt; pf[7] = m.nu;
154+ pf[8] = m.D; pf[9] = m.C; pf[10] = m.G; pf[11] = m.rhoh;
155+ pf[12] = params.Lx / 2; pf[13] = -params.Ly / 2; pf[14] = params.fingerRadius ** 2; pf[15] = mass;
156+ pf[16] = Ix; pf[17] = Iy; pf[18] = Ix + Iy; pf[19] = 0;
157+ this.paramsBuf = device.createBuffer({ size: pf.byteLength, usage: GPUBufferUsage.UNIFORM | GPUBufferUsage.COPY_DST });
158+ device.queue.writeBuffer(this.paramsBuf, 0, pf);
159+
160+ // Per-step parameters, one 32-byte struct per step at a dynamic offset.
161+ this.stepStride = Math.max(32, device.limits.minUniformBufferOffsetAlignment);
162+ this.stepData = new ArrayBuffer(this.stepStride * this.stepsPerBatch);
163+ this.stepBuf = device.createBuffer({ size: this.stepData.byteLength, usage: GPUBufferUsage.UNIFORM | GPUBufferUsage.COPY_DST });
164+
165+ const mk = (code: string, entry: string, layout: GPUPipelineLayout | "auto" = "auto") =>
166+ device.createComputePipeline({
167+ layout,
168+ compute: { module: device.createShaderModule({ code: commonWgsl + code }), entryPoint: entry },
169+ });
170+ const entry = (binding: number, type: GPUBufferBindingType, hasDynamicOffset = false): GPUBindGroupLayoutEntry => ({
171+ binding,
172+ visibility: GPUShaderStage.COMPUTE,
173+ buffer: { type, hasDynamicOffset },
174+ });
175+ const stepLayout = device.createBindGroupLayout({
176+ entries: [
177+ entry(0, "uniform"),
178+ entry(1, "uniform", true),
179+ entry(2, "read-only-storage"),
180+ entry(3, "read-only-storage"),
181+ entry(4, "storage"),
182+ entry(5, "storage"),
183+ entry(6, "storage"),
184+ ],
185+ });
186+ this.stepPipe = mk(stepWgsl, "main", device.createPipelineLayout({ bindGroupLayouts: [stepLayout] }));
187+ const audioLayout = device.createBindGroupLayout({
188+ entries: [entry(0, "uniform"), entry(1, "uniform", true), entry(2, "read-only-storage"), entry(3, "read-only-storage"), entry(4, "storage")],
189+ });
190+ const rigidModule = device.createShaderModule({ code: commonWgsl + rigidWgsl });
191+ const rbLayout = device.createBindGroupLayout({
192+ entries: [
193+ { binding: 0, visibility: GPUShaderStage.COMPUTE, buffer: { type: "uniform" } },
194+ { binding: 1, visibility: GPUShaderStage.COMPUTE, buffer: { type: "storage" } },
195+ { binding: 2, visibility: GPUShaderStage.COMPUTE, buffer: { type: "storage" } },
196+ { binding: 3, visibility: GPUShaderStage.COMPUTE, buffer: { type: "storage" } },
197+ ],
198+ });
199+ const rbPipeLayout = device.createPipelineLayout({ bindGroupLayouts: [rbLayout] });
200+ this.sumsPipe = device.createComputePipeline({ layout: rbPipeLayout, compute: { module: rigidModule, entryPoint: "sums" } });
201+ this.applyPipe = device.createComputePipeline({ layout: rbPipeLayout, compute: { module: rigidModule, entryPoint: "apply" } });
202+ this.audioPipe = mk(audioWgsl, "main", device.createPipelineLayout({ bindGroupLayouts: [audioLayout] }));
203+ this.extractPipe = mk(extractWgsl, "main");
204+
205+ const sp = { buffer: this.stepBuf, size: 32 };
206+ this.stepBG = [0, 1].map((p) =>
207+ device.createBindGroup({
208+ layout: stepLayout,
209+ entries: [
210+ { binding: 0, resource: { buffer: this.paramsBuf } },
211+ { binding: 1, resource: sp },
212+ { binding: 2, resource: { buffer: this.pos[p] } },
213+ { binding: 3, resource: { buffer: this.vel[p] } },
214+ { binding: 4, resource: { buffer: this.pos[1 - p] } },
215+ { binding: 5, resource: { buffer: this.vel[1 - p] } },
216+ { binding: 6, resource: { buffer: this.hist } },
217+ ],
218+ }),
219+ );
220+ this.rbBG = [0, 1].map((p) =>
221+ device.createBindGroup({
222+ layout: rbLayout,
223+ entries: [
224+ { binding: 0, resource: { buffer: this.paramsBuf } },
225+ { binding: 1, resource: { buffer: this.pos[p] } },
226+ { binding: 2, resource: { buffer: this.vel[p] } },
227+ { binding: 3, resource: { buffer: this.statsBuf } },
228+ ],
229+ }),
230+ );
231+ this.audioBG = device.createBindGroup({
232+ layout: audioLayout,
233+ entries: [
234+ { binding: 0, resource: { buffer: this.paramsBuf } },
235+ { binding: 1, resource: sp },
236+ { binding: 2, resource: { buffer: this.hist } },
237+ { binding: 3, resource: { buffer: this.lw } },
238+ { binding: 4, resource: { buffer: this.audioOut } },
239+ ],
240+ });
241+ const extractBG = (p: number, out: GPUBuffer) =>
242+ device.createBindGroup({
243+ layout: this.extractPipe.getBindGroupLayout(0),
244+ entries: [
245+ { binding: 0, resource: { buffer: this.paramsBuf } },
246+ { binding: 1, resource: { buffer: this.pos[p] } },
247+ { binding: 2, resource: { buffer: out } },
248+ ],
249+ });
250+ this.extractBG = [0, 1].map((p) => extractBG(p, this.frameBuf));
251+ this.displayBG = [0, 1].map((p) => extractBG(p, this.displayW));
252+
253+ this.setListener(params.listenerDist, params.listenerAzimuth, params.listenerElevation, params.radiation);
254+ }
255+
256+ /** Recomputes the per-node weights and delays of the listener sum. Takes effect on later samples. */
257+ setListener(dist: number, azimuthDeg: number, elevationDeg: number, radiation: "rayleigh" | "dipole"): void {
258+ const g = this.grid;
259+ const az = (azimuthDeg * Math.PI) / 180;
260+ const el = (elevationDeg * Math.PI) / 180;
261+ const L = [dist * Math.sin(az) * Math.cos(el), dist * Math.sin(el), dist * Math.cos(az) * Math.cos(el)];
262+ const data = new Float32Array(2 * g.n);
263+ let dmin = Infinity;
264+ for (let j = 0; j < g.ny; j++) {
265+ for (let i = 0; i < g.nx; i++) {
266+ const k = j * g.nx + i;
267+ const rx = L[0] - nodeX(g, i), ry = L[1] - nodeY(g, j), rz = L[2];
268+ const r = Math.sqrt(rx * rx + ry * ry + rz * rz);
269+ const dA = g.hx * g.hy * trapWeight(g, i, j);
270+ let wgt = (RHO_AIR * dA) / (2 * Math.PI * r);
271+ if (radiation === "dipole") wgt *= Math.abs(rz) / r;
272+ data[2 * k] = wgt;
273+ data[2 * k + 1] = (r / SPEED_OF_SOUND) * this.params.fs;
274+ dmin = Math.min(dmin, data[2 * k + 1]);
275+ }
276+ }
277+ this.delayClamped = false;
278+ for (let k = 0; k < g.n; k++) {
279+ let d = data[2 * k + 1] - dmin;
280+ if (d > HIST_LEN - 2) {
281+ d = HIST_LEN - 2;
282+ this.delayClamped = true;
283+ }
284+ data[2 * k + 1] = d;
285+ }
286+ this.device.queue.writeBuffer(this.lw, 0, data);
287+ this.params.listenerDist = dist;
288+ this.params.listenerAzimuth = azimuthDeg;
289+ this.params.listenerElevation = elevationDeg;
290+ this.params.radiation = radiation;
291+ }
292+
293+ setDamping(sigma0: number, sigma1: number): void {
294+ this.params.sigma0 = sigma0;
295+ this.params.sigma1 = sigma1;
296+ }
297+
298+ setDuration(duration: number): void {
299+ this.params.duration = duration;
300+ if (this.phase === "done" && this.samples < this.targetSamples()) this.phase = "recording";
301+ }
302+
303+ targetSamples(): number {
304+ return Math.round(this.params.duration * this.params.fs);
305+ }
306+
307+ status(): SimStatus {
308+ const t = this.phase === "push" || this.phase === "hold" ? this.prepTime - this.prepDuration() : this.stepCount * this.dt;
309+ return {
310+ phase: this.phase,
311+ t,
312+ totalSteps: this.totalSteps,
313+ ke: this.ke,
314+ wmax: this.wmax,
315+ samples: this.samples,
316+ frames: this.frames.length,
317+ peakPressure: this.peakPressure,
318+ delayClamped: this.delayClamped,
319+ };
320+ }
321+
322+ private prepDuration(): number {
323+ return this.params.rampTime + this.params.holdMax;
324+ }
325+
326+ /** Enqueues an extraction of the current w field into `displayW` (for the live view). */
327+ encodeDisplay(encoder: GPUCommandEncoder): void {
328+ const pass = encoder.beginComputePass();
329+ pass.setPipeline(this.extractPipe);
330+ pass.setBindGroup(0, this.displayBG[this.cur]);
331+ pass.dispatchWorkgroups(Math.ceil(this.grid.n / 256));
332+ pass.end();
333+ }
334+
335+ /** Listener audio with the radiation post-processing applied [Pa]. */
336+ processedAudio(): Float32Array {
337+ return this.processed.subarray(0, this.samples);
338+ }
339+ private processed: Float32Array = new Float32Array(0);
340+
341+ private appendAudio(raw: Float32Array): void {
342+ const need = this.samples + raw.length;
343+ if (this.audio.length < need) {
344+ const grow = Math.max(need, 2 * this.audio.length, this.params.fs);
345+ const a = new Float32Array(grow);
346+ a.set(this.audio.subarray(0, this.samples));
347+ this.audio = a;
348+ const b = new Float32Array(grow);
349+ b.set(this.processed.subarray(0, this.samples));
350+ this.processed = b;
351+ }
352+ this.audio.set(raw, this.samples);
353+ // Dipole option: first-order high-pass at the baffle cutoff c / (2 sqrt(Lx Ly)).
354+ const dipole = this.params.radiation === "dipole";
355+ const fc = SPEED_OF_SOUND / (2 * Math.sqrt(this.params.Lx * this.params.Ly));
356+ const alpha = 1 / (1 + (2 * Math.PI * fc) / this.params.fs);
357+ for (let k = 0; k < raw.length; k++) {
358+ const x = raw[k];
359+ let y = x;
360+ if (dipole) {
361+ y = alpha * (this.hpState.y + x - this.hpState.x);
362+ }
363+ this.hpState = { y, x };
364+ this.processed[this.samples + k] = y;
365+ this.peakPressure = Math.max(this.peakPressure, Math.abs(y));
366+ }
367+ this.samples = need;
368+ }
369+
370+ /** Runs one batch (one movie frame of time). Resolves when its results have been read back. */
371+ async runBatch(): Promise<void> {
372+ if (this.phase === "done") return;
373+ const p = this.params;
374+ const recording = this.phase === "recording";
375+ const nSteps = this.stepsPerBatch;
376+ const stepsLeft = recording ? (this.targetSamples() - this.samples) * this.stepsPerSample : nSteps;
377+ const steps = Math.min(nSteps, Math.max(0, stepsLeft));
378+ if (steps === 0) {
379+ this.phase = "done";
380+ return;
381+ }
382+
383+ // Fill the per-step parameters.
384+ const u32 = new Uint32Array(this.stepData);
385+ const i32 = new Int32Array(this.stepData);
386+ const f32 = new Float32Array(this.stepData);
387+ const stride = this.stepStride / 4;
388+ const audioStepsInBatch: number[] = [];
389+ for (let s = 0; s < steps; s++) {
390+ const o = s * stride;
391+ const local = this.stepCount + s;
392+ u32[o] = this.totalSteps + s;
393+ let flags = 0, pushW = 0, pushV = 0, sigma0 = p.sigma0, slot = -1, outIdx = -1;
394+ if (!recording) {
395+ const t = this.prepTime + (s + 1) * this.dt; // state after this step
396+ const tau = Math.min(1, t / p.rampTime);
397+ const smooth = tau * tau * (3 - 2 * tau);
398+ const dsmooth = tau < 1 ? (6 * tau * (1 - tau)) / p.rampTime : 0;
399+ flags = 1;
400+ pushW = p.pushAmount * smooth;
401+ pushV = p.pushAmount * dsmooth;
402+ sigma0 = p.settleSigma0;
403+ } else if (local % this.stepsPerSample === 0) {
404+ const sampleIdx = local / this.stepsPerSample;
405+ slot = sampleIdx % HIST_LEN;
406+ outIdx = audioStepsInBatch.length;
407+ audioStepsInBatch.push(s);
408+ }
409+ u32[o + 1] = flags;
410+ i32[o + 2] = slot;
411+ i32[o + 3] = outIdx;
412+ f32[o + 4] = pushW;
413+ f32[o + 5] = pushV;
414+ f32[o + 6] = sigma0;
415+ f32[o + 7] = p.sigma1;
416+ }
417+ this.device.queue.writeBuffer(this.stepBuf, 0, this.stepData, 0, steps * this.stepStride);
418+
419+ // Encode the batch in chunks so that the GPU starts while JS is still encoding.
420+ const CHUNK = 1024;
421+ let enc = this.device.createCommandEncoder();
422+ let pass = enc.beginComputePass();
423+ const wgX = Math.ceil(this.grid.nx / 16), wgY = Math.ceil(this.grid.ny / 16);
424+ const wgN = Math.ceil(this.grid.n / 256);
425+ let cur = this.cur;
426+ let audioPtr = 0;
427+ const rbInterval = recording ? RB_INTERVAL : RB_INTERVAL_HELD;
428+ let stepBound = false;
429+ for (let s = 0; s < steps; s++) {
430+ if (s > 0 && s % CHUNK === 0) {
431+ pass.end();
432+ this.device.queue.submit([enc.finish()]);
433+ enc = this.device.createCommandEncoder();
434+ pass = enc.beginComputePass();
435+ stepBound = false;
436+ }
437+ if (!stepBound) {
438+ pass.setPipeline(this.stepPipe);
439+ stepBound = true;
440+ }
441+ pass.setBindGroup(0, this.stepBG[cur], [s * this.stepStride]);
442+ pass.dispatchWorkgroups(wgX, wgY);
443+ cur = 1 - cur;
444+ const global = this.totalSteps + s + 1;
445+ if (global % rbInterval === 0 || s === steps - 1) {
446+ pass.setPipeline(this.sumsPipe);
447+ pass.setBindGroup(0, this.rbBG[cur]);
448+ pass.dispatchWorkgroups(1);
449+ pass.setPipeline(this.applyPipe);
450+ pass.setBindGroup(0, this.rbBG[cur]);
451+ pass.dispatchWorkgroups(wgN);
452+ stepBound = false;
453+ }
454+ if (audioPtr < audioStepsInBatch.length && audioStepsInBatch[audioPtr] === s) {
455+ pass.setPipeline(this.audioPipe);
456+ pass.setBindGroup(0, this.audioBG, [s * this.stepStride]);
457+ pass.dispatchWorkgroups(1);
458+ audioPtr++;
459+ stepBound = false;
460+ }
461+ }
462+ pass.setPipeline(this.extractPipe);
463+ pass.setBindGroup(0, this.extractBG[cur]);
464+ pass.dispatchWorkgroups(wgN);
465+ pass.end();
466+ const staging = this.staging[this.stagingIdx];
467+ this.stagingIdx = 1 - this.stagingIdx;
468+ const nAudio = audioStepsInBatch.length;
469+ enc.copyBufferToBuffer(this.audioOut, 0, staging, 0, Math.max(4, nAudio * 4));
470+ enc.copyBufferToBuffer(this.frameBuf, 0, staging, this.samplesPerBatch * 4, this.grid.n * 4);
471+ enc.copyBufferToBuffer(this.statsBuf, 0, staging, (this.samplesPerBatch + this.grid.n) * 4, STATS_LEN * 4);
472+ this.device.queue.submit([enc.finish()]);
473+ this.cur = cur;
474+ this.totalSteps += steps;
475+ this.stepCount += steps;
476+ if (!recording) this.prepTime += steps * this.dt;
477+
478+ await staging.mapAsync(GPUMapMode.READ);
479+ const view = new Float32Array(staging.getMappedRange());
480+ const stats = view.subarray(this.samplesPerBatch + this.grid.n, this.samplesPerBatch + this.grid.n + STATS_LEN);
481+ this.ke = stats[12];
482+ this.wmax = stats[13];
483+ if (recording) {
484+ this.appendAudio(view.slice(0, nAudio));
485+ this.frames.push(view.slice(this.samplesPerBatch, this.samplesPerBatch + this.grid.n));
486+ } else {
487+ this.lastPrepFrame = view.slice(this.samplesPerBatch, this.samplesPerBatch + this.grid.n);
488+ }
489+ staging.unmap();
490+
491+ if (!recording) {
492+ if (this.phase === "push" && this.prepTime >= p.rampTime) this.phase = "hold";
493+ const settled = this.phase === "hold" && this.ke < p.settleKE;
494+ if (settled || this.prepTime >= this.prepDuration()) this.release();
495+ } else if (this.samples >= this.targetSamples()) {
496+ this.phase = "done";
497+ }
498+ }
499+ private lastPrepFrame: Float32Array<ArrayBuffer> | null = null;
500+
501+ private release(): void {
502+ this.phase = "recording";
503+ this.stepCount = 0;
504+ this.totalSteps = 0;
505+ // Frame 0 is the held, settled shape.
506+ if (this.lastPrepFrame) this.frames.push(this.lastPrepFrame);
507+ }
508+
509+ /** Reads back the full state (debugging and tests). */
510+ async readState(): Promise<{ pos: Float32Array; vel: Float32Array }> {
511+ const n = this.grid.n;
512+ const stg = this.device.createBuffer({ size: 7 * n * 4, usage: GPUBufferUsage.MAP_READ | GPUBufferUsage.COPY_DST });
513+ const enc = this.device.createCommandEncoder();
514+ enc.copyBufferToBuffer(this.pos[this.cur], 0, stg, 0, 4 * n * 4);
515+ enc.copyBufferToBuffer(this.vel[this.cur], 0, stg, 4 * n * 4, 3 * n * 4);
516+ this.device.queue.submit([enc.finish()]);
517+ await stg.mapAsync(GPUMapMode.READ);
518+ const all = new Float32Array(stg.getMappedRange()).slice();
519+ stg.unmap();
520+ stg.destroy();
521+ // Returned as [u | v | w | wLo] and [du | dv | dw].
522+ return { pos: all.subarray(0, 4 * n), vel: all.subarray(4 * n) };
523+ }
524+
525+ destroy(): void {
526+ for (const b of [...this.pos, ...this.vel, this.hist, this.lw, this.audioOut, this.frameBuf, this.displayW, this.statsBuf, this.paramsBuf, this.stepBuf, ...this.staging]) b.destroy();
527+ }
528+
529+ get fs(): number {
530+ return this.params.fs;
531+ }
532+ get fps(): number {
533+ return this.params.fps;
534+ }
535+ get physParams(): PhysParams {
536+ return this.params;
537+ }
538+}
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; }
tests/cpuRef.test.tsadded+150−0View file
@@ -0,0 +1,150 @@
1+import { test } from "node:test";
2+import assert from "node:assert/strict";
3+import { accelerations, potentialEnergy, maxAngularFrequency, type Fields } from "../src/sim/cpuRef.ts";
4+import { type Grid, type Material, trapWeight } from "../src/params.ts";
5+
6+const mat: Material = { D: 8.66, C: 1.665e8, G: 5.91e7, nu: 0.3, rhoh: 6.241 };
7+
8+function fields(n: number): Fields {
9+ return { u: new Float64Array(n), v: new Float64Array(n), w: new Float64Array(n) };
10+}
11+
12+/** Jacobi eigenvalues of a symmetric matrix (row-major, n x n). */
13+function symEig(A: Float64Array, n: number): number[] {
14+ const a = Float64Array.from(A);
15+ for (let sweep = 0; sweep < 60; sweep++) {
16+ let off = 0;
17+ for (let p = 0; p < n; p++) for (let q = p + 1; q < n; q++) off += a[p * n + q] ** 2;
18+ if (off < 1e-22) break;
19+ for (let p = 0; p < n; p++) {
20+ for (let q = p + 1; q < n; q++) {
21+ const apq = a[p * n + q];
22+ if (Math.abs(apq) < 1e-300) continue;
23+ const theta = (a[q * n + q] - a[p * n + p]) / (2 * apq);
24+ const t = Math.sign(theta || 1) / (Math.abs(theta) + Math.sqrt(theta * theta + 1));
25+ const c = 1 / Math.sqrt(t * t + 1), s = t * c;
26+ for (let k = 0; k < n; k++) {
27+ const akp = a[k * n + p], akq = a[k * n + q];
28+ a[k * n + p] = c * akp - s * akq;
29+ a[k * n + q] = s * akp + c * akq;
30+ }
31+ for (let k = 0; k < n; k++) {
32+ const apk = a[p * n + k], aqk = a[q * n + k];
33+ a[p * n + k] = c * apk - s * aqk;
34+ a[q * n + k] = s * apk + c * aqk;
35+ }
36+ }
37+ }
38+ }
39+ const ev: number[] = [];
40+ for (let k = 0; k < n; k++) ev.push(a[k * n + k]);
41+ return ev.sort((x, y) => x - y);
42+}
43+
44+/** Mass-symmetrized linear stiffness of one block ("w" or "uv"). */
45+function stiffness(g: Grid, block: "w" | "uv"): { K: Float64Array; n: number } {
46+ const N = g.n;
47+ const n = block === "w" ? N : 2 * N;
48+ const K = new Float64Array(n * n);
49+ const x = fields(N), y = fields(N);
50+ const eps = 1e-9;
51+ const mass = (k: number) => mat.rhoh * g.hx * g.hy * trapWeight(g, k % g.nx, Math.floor(k / g.nx));
52+ for (let col = 0; col < n; col++) {
53+ x.u.fill(0); x.v.fill(0); x.w.fill(0);
54+ if (block === "w") x.w[col] = eps;
55+ else if (col < N) x.u[col] = eps;
56+ else x.v[col - N] = eps;
57+ accelerations(g, mat, x, y);
58+ const mc = mass(col % N);
59+ for (let row = 0; row < n; row++) {
60+ const acc = block === "w" ? y.w[row] : row < N ? y.u[row] : y.v[row - N];
61+ const mr = mass(row % N);
62+ // K_sym = M^{1/2} (M^{-1} K) M^{-1/2}
63+ K[row * n + col] = (-acc / eps) * Math.sqrt(mr / mc);
64+ }
65+ }
66+ return { K, n };
67+}
68+
69+test("bending stiffness is symmetric and has exactly three rigid modes", () => {
70+ const nx = 9, ny = 7;
71+ const g: Grid = { nx, ny, hx: 0.3 / (nx - 1), hy: 0.2 / (ny - 1), n: nx * ny };
72+ const { K, n } = stiffness(g, "w");
73+ let asym = 0, scale = 0;
74+ for (let r = 0; r < n; r++) for (let c = 0; c < n; c++) {
75+ asym = Math.max(asym, Math.abs(K[r * n + c] - K[c * n + r]));
76+ scale = Math.max(scale, Math.abs(K[r * n + c]));
77+ }
78+ assert.ok(asym < 1e-6 * scale, `asymmetry ${asym} vs scale ${scale}`);
79+ const ev = symEig(K, n);
80+ assert.ok(Math.abs(ev[2]) < 1e-6 * ev[n - 1], "third rigid mode");
81+ assert.ok(ev[3] > 1e-4 * ev[n - 1], "fourth mode is not rigid");
82+});
83+
84+test("membrane stiffness is symmetric and has exactly three rigid modes", () => {
85+ const nx = 7, ny = 6;
86+ const g: Grid = { nx, ny, hx: 0.3 / (nx - 1), hy: 0.2 / (ny - 1), n: nx * ny };
87+ const { K, n } = stiffness(g, "uv");
88+ let asym = 0, scale = 0;
89+ for (let r = 0; r < n; r++) for (let c = 0; c < n; c++) {
90+ asym = Math.max(asym, Math.abs(K[r * n + c] - K[c * n + r]));
91+ scale = Math.max(scale, Math.abs(K[r * n + c]));
92+ }
93+ assert.ok(asym < 1e-6 * scale, `asymmetry ${asym}`);
94+ const ev = symEig(K, n);
95+ assert.ok(Math.abs(ev[2]) < 1e-6 * ev[n - 1], "third rigid mode");
96+ assert.ok(ev[3] > 1e-4 * ev[n - 1], "fourth mode is not rigid");
97+});
98+
99+test("free square plate frequencies match the Leissa values", () => {
100+ const nx = 21, ny = 21, a = 0.3;
101+ const g: Grid = { nx, ny, hx: a / (nx - 1), hy: a / (ny - 1), n: nx * ny };
102+ const { K, n } = stiffness(g, "w");
103+ const ev = symEig(K, n);
104+ const lam2 = ev.slice(3, 8).map((e) => Math.sqrt(e) * a * a * Math.sqrt(mat.rhoh / mat.D));
105+ // Completely free square plate, nu = 0.3 (Leissa / Gorman): first five nonzero modes.
106+ const ref = [13.468, 19.596, 24.27, 34.801, 34.801];
107+ for (let k = 0; k < ref.length; k++) {
108+ const relErr = Math.abs(lam2[k] - ref[k]) / ref[k];
109+ assert.ok(relErr < 0.03, `mode ${k}: got ${lam2[k].toFixed(3)}, expected ${ref[k]} (err ${(relErr * 100).toFixed(2)}%)`);
110+ }
111+});
112+
113+test("forces are the negative gradient of the discrete energy (nonlinear state)", () => {
114+ const nx = 8, ny = 7;
115+ const g: Grid = { nx, ny, hx: 0.3 / (nx - 1), hy: 0.2 / (ny - 1), n: nx * ny };
116+ const f = fields(g.n), acc = fields(g.n);
117+ let seed = 7;
118+ const rnd = () => ((seed = (seed * 1664525 + 1013904223) >>> 0) / 4294967296 - 0.5);
119+ for (let k = 0; k < g.n; k++) {
120+ f.w[k] = 0.02 * rnd(); // ~ 25 thicknesses, strongly nonlinear
121+ f.u[k] = 1e-4 * rnd();
122+ f.v[k] = 1e-4 * rnd();
123+ }
124+ accelerations(g, mat, f, acc);
125+ const delta = 1e-7;
126+ let worst = 0;
127+ for (const [arr, accArr] of [[f.w, acc.w], [f.u, acc.u], [f.v, acc.v]] as const) {
128+ for (let k = 0; k < g.n; k += 5) {
129+ const mass = mat.rhoh * g.hx * g.hy * trapWeight(g, k % nx, Math.floor(k / nx));
130+ const orig = arr[k];
131+ arr[k] = orig + delta;
132+ const Ep = potentialEnergy(g, mat, f);
133+ arr[k] = orig - delta;
134+ const Em = potentialEnergy(g, mat, f);
135+ arr[k] = orig;
136+ const Ffd = -(Ep - Em) / (2 * delta);
137+ const Fan = accArr[k] * mass;
138+ worst = Math.max(worst, Math.abs(Ffd - Fan) / (Math.abs(Fan) + 1e-3));
139+ }
140+ }
141+ assert.ok(worst < 1e-4, `gradient mismatch ${worst}`);
142+});
143+
144+test("maximum frequency estimate is close to the in-plane CFL estimate", () => {
145+ const g: Grid = { nx: 50, ny: 50, hx: 0.004, hy: 0.004, n: 2500 };
146+ const wmax = maxAngularFrequency(g, mat);
147+ const cp = Math.sqrt(mat.C / mat.rhoh);
148+ const cfl = (2 * cp * Math.SQRT2) / 0.004; // 2 c sqrt(1/hx^2 + 1/hy^2)
149+ assert.ok(wmax > 0.5 * cfl && wmax < 1.5 * cfl, `wmax ${wmax} vs cfl ${cfl}`);
150+});
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"]
17+ },
18+ "include": ["src"]
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+});