Drumstick by default, hit sequences, 170 GPa steel
The drumstick strike is the default excitation. A Hits selector turns the
single hit into a series (same speed, soft to loud, loud to soft) with a
hit count and interval; at each scheduled step the shader re-arms the
stick at the current contact coordinate with the scheduled speed. The
Excitation panel reports the schedule and the ring left after the last
hit. Default Young's modulus 170 GPa.
Also: in-plane mass factor (scaled in-plane inertia, lengthens the time
step), one audio dispatch per batch, tiled step kernel with shared-memory
stencil ingredients, and a self-test case for the mass factor.
14 changed files+434−152
README.mdmodified+4−4View file
@@ -17,7 +17,7 @@ Requires Node 20 or later and a browser with WebGPU (Chrome or Edge 113+; Firefo
1717
1818 ## Using the simulator
1919
20-Press Start. With the push excitation the sheet is pushed at the bottom-right corner over the ramp time, held with heavy artificial damping until it settles, and released by ramping the finger force to zero over the release time; the recording of audio and movie frames starts at the release. With the drumstick excitation a free mass strikes the sheet at a configurable point and speed through a Hertz contact, and the recording starts at the impact. The view shows the transverse displacement as a color map (2D) or a lit height field (3D; drag to orbit, wheel to zoom), the spectrogram shows the listener signal, and the audio can be played or downloaded as a 16-bit WAV at any time.
20+Press Start. With the drumstick excitation (the default) a free mass strikes the sheet at a configurable point and speed through a Hertz contact, and the recording of audio and movie frames starts at the impact. The Hits selector turns the single hit into a series: a fixed number of hits at a fixed interval, all at the same speed or rising or falling between the softest and the loudest speed; the Excitation panel reports the schedule and how much ring the recorded duration leaves after the last hit. With the push excitation the sheet is pushed at the bottom-right corner over the ramp time, held with heavy artificial damping until it settles (its kinetic energy drops below the threshold or stops falling), and released by ramping the finger force to zero over the release time; the recording starts at the release. The view shows the transverse displacement as a color map (2D) or a lit height field (3D; drag to orbit, wheel to zoom), the spectrogram shows the listener signal, and the audio can be played or downloaded as a 16-bit WAV at any time.
2121
2222 - The sheet is a rectangle by default; the Outline selector makes it the ellipse inscribed in the given width and height. Width and height are entered in inches by default, or in millimeters; everything else is metric. Internally all values are SI.
2323 - The sheet has a stress-free bow (default 3 mm across the width; set the bow shape to Flat for the ideal plate). A negative push amount pushes against the bow.
@@ -26,16 +26,16 @@ Press Start. With the push excitation the sheet is pushed at the bottom-right co
2626 - Pause halts the run temporarily; Stop ends it, keeps the audio and movie, and unlocks the parameter fields. To extend a finished run, raise the duration and press Continue; the state is kept on the GPU.
2727 - Playback level: Absolute maps 1 Pa to full scale, shifted by the gain slider, so different runs can be compared by ear (the peak level in dBFS is shown, with a clipping warning); Loudness and Peak normalize each run instead. The WAV download uses the same setting.
2828 - After the simulation completes, drag the time slider to explore the movie; Play starts the audio from the slider position and moves the movie and the spectrogram cursor with it. Live returns to the running view.
29-- Any parameter can be preset in the URL in SI units, for example `?spacing=0.006&duration=2&autostart=1` or `?excitation=strike&shape=ellipse`.
29+- Any parameter can be preset in the URL in SI units, for example `?spacing=0.006&duration=2&autostart=1` or `?excitation=push&shape=ellipse` or `?hitSequence=crescendo&hitCount=6&hitInterval=0.5&duration=6`.
3030
31-The default grid spacing of 8 mm (39 by 58 nodes) keeps the simulation interactive but lowers the bending partials by 3 percent near 1 kHz and by 10 percent near 3.6 kHz through grid dispersion; the Sheet panel reports these two frequencies for the chosen spacing. The stable time step is set by the in-plane wave speed, about 1.3 microseconds at 8 mm, so 4 s of sound take about 3 million steps. On an Intel Iris Xe integrated GPU this runs at about 33k steps per second, or roughly 1.5 minutes for 4 s of sound; a discrete GPU is considerably faster. Halving the spacing to 4 mm costs about eight times as much (about 5 minutes for 4 s on the same GPU) and moves the 3 percent point to about 4 kHz.
31+The default grid spacing of 8 mm (39 by 58 nodes) keeps the simulation interactive but lowers the bending partials by 3 percent near 1 kHz and by 10 percent near 3.6 kHz through grid dispersion; the Sheet panel reports these two frequencies for the chosen spacing. The stable time step is set by the in-plane wave speed, about 1.4 microseconds at 8 mm, so 4 s of sound take about 2.9 million steps. On an Intel Iris Xe integrated GPU this runs at about 57k steps per second, or roughly 50 s for 4 s of sound plus a few seconds for the push and hold; a discrete GPU is considerably faster. On grids of this size the cost of a step is the latency of one workgroup rather than the arithmetic, so a finer grid costs less than the node count suggests: halving the spacing to 4 mm costs about three times as much (about 2.5 minutes for 4 s on the same GPU) and moves the 3 percent point to about 4 kHz. The Numerics panel also offers an in-plane mass factor that lengthens the time step by its square root at the expense of lowering the in-plane modes of the sheet; see THEORY.md, Section 5.
3232
3333 ## Tests
3434
3535 `npm test` runs two groups of tests with Node's built-in runner.
3636
3737 - `tests/cpuRef.test.ts` checks the float64 reference of the spatial discretization without a GPU: symmetry of the stiffness with exactly three rigid modes per block, forces as the exact gradient of the discrete energy (flat and bowed, rectangle and ellipse), dissipativity of the loss terms, and convergence of the free square and circular plate frequencies to the reference values of Gorman and Leissa.
38-- `tests/gpu.test.ts` serves the app with Vite and drives a headless Chromium through Playwright. It compares one step of the WGSL kernel with the float64 reference on flat, bowed, and elliptical sheets (`?selftest=1` in the browser shows the same numbers), and runs a short push and a short strike to completion. It needs Playwright's Chromium (`npx playwright install chromium`) and a machine whose Chromium has WebGPU; `npm run test:cpu` runs only the first group.
38+- `tests/gpu.test.ts` serves the app with Vite and drives a headless Chromium through Playwright. It compares one step of the WGSL kernel with the float64 reference on flat, bowed, and elliptical sheets (`?selftest=1` in the browser shows the same numbers), and runs a short push, a short strike, and a short series of strikes to completion. It needs Playwright's Chromium (`npx playwright install chromium`) and a machine whose Chromium has WebGPU; `npm run test:cpu` runs only the first group.
3939
4040 ## Layout
4141
THEORY.mdmodified+8−6View file
@@ -38,13 +38,13 @@ The whole edge is free: the bending moment and the Kirchhoff shear vanish, and s
3838
3939 The thread through the hole near one corner is not modeled as a constraint. At audio frequencies it is a very soft spring compared with the local stiffness of the plate, so it does not pin the hole; it only fixes the average position of the plate. We therefore treat the plate as free, and we make every external load self-equilibrated: the reaction to the finger or the stick is distributed over the plate as the inertial load of a rigid body, which is the inertia-relief equilibrium of a free body. The sheet then gains no net momentum or angular momentum, which stands in for the thread taking them up. This is a modeling choice: a sheet hanging freely from a thread and pushed 10 mm at a corner would mostly swing as a pendulum and bend by only a millimeter or two, and a struck sheet would recoil and swing. Gravity is dropped: the in-plane stress it induces is about $\rho g L \approx 35$ kPa, a strain of $2\times10^{-7}$.
4040
41-Inertia relief takes two forms here. During a push the finger prescribes the displacement of a small patch, and the reaction is not known in advance, so the six rigid-body modes (three translations, three rotations) are removed from the displacement and velocity fields by a mass-weighted projection every four time steps while the sheet is held; the push therefore produces pure bending of the prescribed amount. The drumstick force is known at every step, so its reaction is applied directly as the rigid acceleration field that balances the force and its moments about the center (Section 3). Once the sheet rings freely, no load remains that could move it as a whole, apart from rounding and the small non-invariance of the von Karman energy under finite rotation, and the projection is optional there and off by default. Note that a rigid rotation is not exactly a null mode of the von Karman energy: the strains contain $\tfrac12 w_x^2$, so a rotation $\theta$ about an in-plane axis produces a spurious strain $\theta^2/2$. This is the reason the loads have to be self-equilibrated, since a slow drift in rotation would otherwise stress the sheet.
41+Inertia relief takes two forms here. During a push the finger prescribes the displacement of a small patch, and the reaction is not known in advance, so the six rigid-body modes (three translations, three rotations) are removed from the displacement and velocity fields by a mass-weighted projection every 16 time steps (about 20 microseconds) while the sheet is held; the push therefore produces pure bending of the prescribed amount. The drumstick force is known at every step, so its reaction is applied directly as the rigid acceleration field that balances the force and its moments about the center (Section 3). Once the sheet rings freely, no load remains that could move it as a whole, apart from rounding and the small non-invariance of the von Karman energy under finite rotation, and the projection is optional there and off by default. Note that a rigid rotation is not exactly a null mode of the von Karman energy: the strains contain $\tfrac12 w_x^2$, so a rotation $\theta$ about an in-plane axis produces a spurious strain $\theta^2/2$. This is the reason the loads have to be self-equilibrated, since a slow drift in rotation would otherwise stress the sheet.
4242
4343 ## 3. Excitation
4444
45-**Push and release.** The simulation has three phases. During the push (default 0.1 s) the transverse displacement of the nodes within the finger radius follows a smooth ramp to the push amount, and heavy artificial damping ($\sigma_0 = 300\ \mathrm{s^{-1}}$) removes the waves that the ramp launches. During the hold the constraint stays in place until the non-rigid kinetic energy drops below a threshold or a maximum settle time elapses. At the release the clock is set to zero, the damping returns to its physical value, and the recording of audio and movie frames begins. The finger does not let go instantaneously: the reaction force that held each node of the patch is recorded at the end of the hold and then applied as an external force that falls smoothly to zero over the release time (20 ms by default). This matters for the sound. An instantaneous release is a step in the finger force, which excites all modes at once, and since the radiated pressure is proportional to the plate acceleration, the result is a nearly flat pressure spectrum, a cymbal-like click. A release over 20 ms confines the excitation to roughly the first hundred hertz, as a fingertip does. The rigid-body components of the recorded reaction are removed before it is applied, so that the finger force and the distributed holding reaction relax together. The release time can be set to zero to recover the instantaneous case.
45+**Push and release.** The simulation has three phases. During the push (default 0.1 s) the transverse displacement of the nodes within the finger radius follows a smooth ramp to the push amount, and heavy artificial damping ($\sigma_0 = 300\ \mathrm{s^{-1}}$) removes the waves that the ramp launches. During the hold the constraint stays in place until the sheet has settled or a maximum settle time elapses. The sheet counts as settled when its non-rigid kinetic energy drops below a threshold, or when that energy has stopped falling: in single precision the rounding of the force evaluation leaves a floor of a few times $10^{-8}$ J that can lie above any fixed threshold, so the hold also ends once the energy has not halved over the last 50 ms (or over three settle time constants, whichever is longer). With the default settle damping the release comes about 0.14 s into the hold. At the release the clock is set to zero, the damping returns to its physical value, and the recording of audio and movie frames begins. The finger does not let go instantaneously: the reaction force that held each node of the patch is recorded at the end of the hold and then applied as an external force that falls smoothly to zero over the release time (20 ms by default). This matters for the sound. An instantaneous release is a step in the finger force, which excites all modes at once, and since the radiated pressure is proportional to the plate acceleration, the result is a nearly flat pressure spectrum, a cymbal-like click. A release over 20 ms confines the excitation to roughly the first hundred hertz, as a fingertip does. The rigid-body components of the recorded reaction are removed before it is applied, so that the finger force and the distributed holding reaction relax together. The release time can be set to zero to recover the instantaneous case.
4646
47-**Drumstick strike.** The stick tip is a free mass $M$ (default 30 g, the effective mass of a stick in a wrist stroke rather than its full mass) arriving with speed $v_0$ (default 3 m/s) at a point of the sheet. The contact patch consists of the sheet nodes within the tip radius (5 mm) plus the nearest node. Its contact coordinate is the mass-weighted mean transverse displacement $w_c$ of these nodes, the penetration is $\delta = w_c - x$ with $x$ the tip position, and the force is the Hertz law $F = K\,\delta^{\alpha}$ for $\delta > 0$ and zero otherwise, with $\alpha = 1.5$ for a rounded tip and $K = 5\times10^{8}\ \mathrm{N/m^{1.5}}$ for wood on steel. The force is spread over the patch in proportion to the node masses. Because that distribution is exactly conjugate to the mean displacement, the contact has a potential energy and the coupled system is Hamiltonian apart from the losses. At the same time the force is balanced over the whole sheet by the inertia-relief field $a_k = F\,(1/M_s + y_k y_c/I_x + x_k x_c/I_y)$, where $M_s$ is the mass of the sheet, $I_x$ and $I_y$ its moments of inertia about the in-plane axes, and $(x_c, y_c)$ the mass-weighted center of the patch, so that the total force and the two moments of the load vanish (Section 2). The field uses the force of the previous step, a lag of one microsecond whose net impulse over the contact is zero. The stick obeys $M\ddot x = F$ and is stepped together with the plate; it may bounce and strike again if the sheet comes back to meet it. The recording starts at the impact, with the sheet at rest in its natural shape. Note that a thin sheet is a soft target: its driving-point impedance is $8\sqrt{D\rho h} \approx 59\ \mathrm{N\,s/m}$, so the contact lasts well under a millisecond and the peak force is of the order of a hundred newtons. The excitation is therefore broadband, unlike the push, and the sound is a short bright attack followed by the ring.
47+**Drumstick strike.** The stick tip is a free mass $M$ (default 30 g, the effective mass of a stick in a wrist stroke rather than its full mass) arriving with speed $v_0$ (default 3 m/s) at a point of the sheet. The contact patch consists of the sheet nodes within the tip radius (5 mm) plus the nearest node. Its contact coordinate is the mass-weighted mean transverse displacement $w_c$ of these nodes, the penetration is $\delta = w_c - x$ with $x$ the tip position, and the force is the Hertz law $F = K\,\delta^{\alpha}$ for $\delta > 0$ and zero otherwise, with $\alpha = 1.5$ for a rounded tip and $K = 5\times10^{8}\ \mathrm{N/m^{1.5}}$ for wood on steel. The force is spread over the patch in proportion to the node masses. Because that distribution is exactly conjugate to the mean displacement, the contact has a potential energy and the coupled system is Hamiltonian apart from the losses. At the same time the force is balanced over the whole sheet by the inertia-relief field $a_k = F\,(1/M_s + y_k y_c/I_x + x_k x_c/I_y)$, where $M_s$ is the mass of the sheet, $I_x$ and $I_y$ its moments of inertia about the in-plane axes, and $(x_c, y_c)$ the mass-weighted center of the patch, so that the total force and the two moments of the load vanish (Section 2). The field uses the force of the previous step, a lag of one microsecond whose net impulse over the contact is zero. The stick obeys $M\ddot x = F$ and is stepped together with the plate; it may bounce and strike again if the sheet comes back to meet it. The recording starts at the impact, with the sheet at rest in its natural shape. A series of hits is scheduled at fixed intervals with speeds that are constant, rising, or falling from hit to hit. At each scheduled instant the stick is replaced: the tip is put at the current contact coordinate $w_c$ of the moving sheet, wherever that is, with the scheduled speed toward the sheet and zero force, and the rebounding stick of the previous hit is discarded. Each hit therefore starts from touching contact exactly as the first does, and the sheet keeps whatever motion it has from the earlier hits. Note that a thin sheet is a soft target: its driving-point impedance is $8\sqrt{D\rho h} \approx 59\ \mathrm{N\,s/m}$, so the contact lasts well under a millisecond and the peak force is of the order of a hundred newtons. The excitation is therefore broadband, unlike the push, and the sound is a short bright attack followed by the ring.
4848
4949 ## 4. Spatial discretization
5050
@@ -64,11 +64,13 @@ The plate is covered by a regular grid of $n_x \times n_y$ nodes with spacings $
6464
6565 ## 5. Time integration
6666
67-The system is advanced with the symplectic Euler scheme, $v^{n+1} = v^n + \Delta t\, a(x^n, v^n)$ and $x^{n+1} = x^n + \Delta t\, v^{n+1}$, in single precision on the GPU. With the velocities interpreted at half steps this is the leapfrog (Stormer-Verlet) scheme, second order and symplectic for the conservative part; the losses are evaluated explicitly at the old velocity. The stability limit is $\Delta t < 2/\omega_{\max}$, where $\omega_{\max}$ is the largest frequency of the linearized discrete system. We estimate $\omega_{\max}$ at start-up by power iteration on a small grid with the same spacing and use 85 percent of the limit. The in-plane waves, with speed $c_p = \sqrt{E/\rho(1-\nu^2)} \approx 5200$ m/s, set this limit: at 8 mm spacing $\Delta t \approx 1.3\ \mu$s, roughly six times smaller than bending alone would require, and 4 s of sound take about 3 million steps; at 4 mm the step halves and the node count quadruples. The step is rounded down so that an integer number of steps (16 at 8 mm and 48 kHz) makes one audio sample.
67+The system is advanced with the symplectic Euler scheme, $v^{n+1} = v^n + \Delta t\, a(x^n, v^n)$ and $x^{n+1} = x^n + \Delta t\, v^{n+1}$, in single precision on the GPU. With the velocities interpreted at half steps this is the leapfrog (Stormer-Verlet) scheme, second order and symplectic for the conservative part; the losses are evaluated explicitly at the old velocity. The stability limit is $\Delta t < 2/\omega_{\max}$, where $\omega_{\max}$ is the largest frequency of the linearized discrete system. We estimate $\omega_{\max}$ at start-up by power iteration on a small grid with the same spacing and use 90 percent of the limit. The in-plane waves, with speed $c_p = \sqrt{E/\rho(1-\nu^2)} \approx 5200$ m/s, set this limit: at 8 mm spacing $\Delta t \approx 1.3\ \mu$s, roughly six times smaller than bending alone would require, and 4 s of sound take about 3 million steps; at 4 mm the step halves and the node count quadruples. The step is rounded down so that an integer number of steps (15 at 8 mm and 48 kHz) makes one audio sample.
68+
69+**In-plane mass scaling.** Since the in-plane waves set the step and the in-plane inertia is a small effect (the usual Airy stress function formulation drops it altogether), the interface offers a factor $s \geq 1$ by which the in-plane inertia is scaled; the default is 1. The in-plane wave speed becomes $c_p/\sqrt{s}$ and the stable step grows by up to $\sqrt{s}$, until bending sets the limit at about $s = 30$ on the 8 mm grid. The trade-off is that the in-plane modes of the sheet, which start near 8 kHz for the default sheet, move down by the same factor and, once the sheet is curved, are coupled to the transverse motion and to the sound. Note that the scaling moves these modes in the opposite direction from the quasi-static limit, which places them at infinite frequency. A factor of 4 halves the run time and keeps the in-plane modes above 4 kHz; it is meant for quick previews at a finer spacing rather than for the final sound. The kinetic energy reported by the interface includes the factor, so the settle criterion is unaffected.
6870
6971 Single precision needs one precaution. With $w$ of a few millimeters and $\Delta t \approx 1\ \mu$s, a velocity below about $10^{-3}$ m/s produces an increment smaller than the resolution of $w$, and the rounding would act as a broadband noise source. The update of $w$ is therefore compensated: the rounding error of each addition is carried in a separate array and added back at the next step (Kahan summation). The self-test of Section 4 also checks this update.
7072
71-The state is stored twice and each step reads one copy and writes the other. A step is one compute dispatch of 16 by 16 workgroups; a per-step uniform, addressed by a dynamic offset, carries the prescribed patch displacement, the damping, and the audio bookkeeping, so that a whole movie frame of steps is encoded in one pass. The mass-weighted sums over the sheet (rigid-body components, kinetic energy, peak displacement) run in one workgroup every four steps while the sheet is held and once per frame otherwise, or every 32 steps when the optional projection of Section 2 is on.
73+The state is stored twice and each step reads one copy and writes the other. A step is one compute dispatch in which each workgroup of 256 threads advances a 12 by 12 tile of nodes: the threads load the tile with a two-node halo into shared memory, compute the membrane forces of the cells and the bending moments at the nodes once each, and then evaluate the update of the tile's nodes from shared memory. On grids of a few thousand nodes the cost of a step is the latency of one such workgroup, about 17 microseconds on an integrated GPU, rather than the arithmetic, so the shared phases, which avoid recomputing the stencil ingredients at every node, roughly halve the step time; the step time barely grows until the grid exceeds several thousand nodes. A per-step uniform, addressed by a dynamic offset, carries the prescribed patch displacement, the damping, and the audio bookkeeping, so that a whole movie frame of steps is encoded in one pass, and the following frame is encoded while the GPU works on the current one. The mass-weighted sums over the sheet (rigid-body components, kinetic energy, peak displacement) run in one workgroup every 16 steps while the sheet is held and once per frame otherwise, or every 32 steps when the optional projection of Section 2 is on.
7274
7375 ## 6. Radiation to the listener
7476
@@ -78,7 +80,7 @@ $$
7880 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),
7981 $$
8082
81-where the sum runs over the nodes, $r_k$ is the distance from node $k$ to the listener, and $c = 343$ m/s. The transverse velocity of every node is recorded in a ring buffer at the audio rate, and each output sample is a delayed, weighted sum over the plate, evaluated on the GPU with linear interpolation between samples; the acceleration is formed as the difference between successive delayed velocities, which is the average acceleration over one sample interval. This matters in single precision: the force evaluation carries a broadband rounding error of the order of $1\ \mathrm{m/s^2}$ per node, and sampling the instantaneous acceleration would put that noise straight into the audio, whereas the averaged form suppresses it by the square root of the number of steps per sample. The common delay is discarded; only the spread of delays across the plate (about 5 samples for a listener 1 m in front) is kept.
83+where the sum runs over the nodes, $r_k$ is the distance from node $k$ to the listener, and $c = 343$ m/s. The transverse velocity of every node is recorded in a ring buffer at the audio rate, and each output sample is a delayed, weighted sum over the plate, evaluated on the GPU once per movie frame for all the frame's samples, with linear interpolation between samples; the acceleration is formed as the difference between successive delayed velocities, which is the average acceleration over one sample interval. This matters in single precision: the force evaluation carries a broadband rounding error of the order of $1\ \mathrm{m/s^2}$ per node, and sampling the instantaneous acceleration would put that noise straight into the audio, whereas the averaged form suppresses it by the square root of the number of steps per sample. The common delay is discarded; only the spread of delays across the plate (about 5 samples for a listener 1 m in front) is kept.
8284
8385 The plate is unbaffled in reality, so the sound from the two faces partly cancels at low frequencies, and the Rayleigh integral overestimates the low-frequency level. An accurate treatment would require solving the exterior acoustic problem, which is outside the scope of this simulator. The pressure is reported in pascals and the peak sound pressure level is displayed. Playback and the WAV file use, by default, a fixed mapping from pressure to full scale (1 Pa at 0 dBFS, shifted by a user gain in decibels), so that runs with different excitations can be compared by ear; a soft clipper bounds the output and the resulting peak level is shown. Two normalized modes are also offered, to the peak or to the loudness of the sustained sound.
8486
src/main.tsmodified+2−0View file
@@ -63,6 +63,7 @@ function buildControls(): void {
6363 duration = v;
6464 params.duration = v;
6565 sim?.setDuration(v);
66+ form.updateNotes();
6667 updateButtons();
6768 }
6869 });
@@ -397,6 +398,7 @@ function updateInfo(st: SimStatus): void {
397398 ["Peak pressure", st.peakPressure > 0 ? `${st.peakPressure.toFixed(3)} Pa (${spl.toFixed(1)} dB SPL)` : "–"],
398399 ];
399400 if (st.stick) {
401+ if (st.hits && st.hits.total > 1) rows.push(["Hits", `${st.hits.done} of ${st.hits.total}`]);
400402 rows.push(["Stick", `${(st.stick.x * 1e3).toFixed(2)} mm, ${st.stick.v.toFixed(2)} m/s`]);
401403 rows.push(["Contact force", `${st.stick.force.toFixed(1)} N (max ${st.stick.maxForce.toFixed(1)} N)`]);
402404 }
src/params.tsmodified+52−5View file
@@ -5,6 +5,7 @@
55 export type PlateShape = "rectangle" | "ellipse";
66 export type BowShape = "none" | "width" | "height" | "dish";
77 export type Excitation = "push" | "strike";
8+export type HitSequence = "single" | "steady" | "crescendo" | "decrescendo";
89
910 export interface PhysParams {
1011 // Sheet
@@ -34,7 +35,11 @@ export interface PhysParams {
3435 settleKE: number; // non-rigid kinetic energy below which the sheet counts as settled [J]
3536 strikeFracX: number; // strike point as a fraction of the width from the left edge
3637 strikeFracY: number; // strike point as a fraction of the height from the bottom edge
37- strikeVelocity: number; // stick tip speed at impact [m/s]
38+ strikeVelocity: number; // stick tip speed at impact [m/s]; the loudest hit of a sequence
39+ hitSequence: HitSequence; // one hit, or a series of hits every hitInterval seconds
40+ hitCount: number; // number of hits in a series
41+ hitInterval: number; // time between the hits of a series [s]
42+ strikeVelocityMin: number; // softest stick speed of a crescendo or decrescendo [m/s]
3843 stickMass: number; // effective mass of the stick [kg]
3944 contactStiffness: number; // Hertz contact stiffness K [N/m^alpha]
4045 contactExponent: number; // Hertz exponent alpha (1.5 for a rounded tip)
@@ -45,6 +50,7 @@ export interface PhysParams {
4550 listenerElevation: number; // [deg]
4651 // Numerics and output
4752 spacing: number; // requested grid spacing [m]
53+ inPlaneMass: number; // factor by which the in-plane inertia is scaled; 1 is physical, larger lengthens the time step
4854 rigidProjection: boolean; // remove the six rigid-body modes every RB_INTERVAL steps while ringing
4955 fs: number; // audio sample rate [Hz]
5056 fps: number; // movie frame rate, also the batch rate [Hz]
@@ -56,7 +62,7 @@ export const defaultParams: PhysParams = {
5662 Lx: 12 * 0.0254,
5763 Ly: 18 * 0.0254,
5864 h: 0.79e-3,
59- E: 193e9,
65+ E: 170e9,
6066 nu: 0.29,
6167 rho: 7900,
6268 holeInset: 0.012,
@@ -65,7 +71,7 @@ export const defaultParams: PhysParams = {
6571 sigma0: 1.15,
6672 sigma1: 5e-4,
6773 sigmaIn: 20,
68- excitation: "push",
74+ excitation: "strike",
6975 pushAmount: 0.01,
7076 fingerRadius: 0.01,
7177 rampTime: 0.1,
@@ -76,6 +82,10 @@ export const defaultParams: PhysParams = {
7682 strikeFracX: 0.7,
7783 strikeFracY: 0.35,
7884 strikeVelocity: 3,
85+ hitSequence: "single",
86+ hitCount: 6,
87+ hitInterval: 0.5,
88+ strikeVelocityMin: 0.5,
7989 stickMass: 0.03,
8090 contactStiffness: 5e8,
8191 contactExponent: 1.5,
@@ -84,6 +94,7 @@ export const defaultParams: PhysParams = {
8494 listenerAzimuth: 0,
8595 listenerElevation: 0,
8696 spacing: 0.008,
97+ inPlaneMass: 1,
8798 rigidProjection: false,
8899 fs: 48000,
89100 fps: 120,
@@ -92,9 +103,17 @@ export const defaultParams: PhysParams = {
92103
93104 export const SPEED_OF_SOUND = 343; // [m/s]
94105 export const RHO_AIR = 1.2; // [kg/m^3]
95-export const HIST_LEN = 256; // ring buffer length of the retarded-time listener sum [samples]
106+export const CFL = 0.9; // fraction of the stability limit used for the time step
107+// The ring buffer of the retarded-time listener sum holds one batch of samples
108+// plus this margin, which bounds the spread of the delays across the sheet.
109+export const HIST_MARGIN = 128; // [samples]
96110 export const RB_INTERVAL = 32; // steps between rigid-body projections while ringing
97-export const RB_INTERVAL_HELD = 4; // while the finger constraint is active
111+export const RB_INTERVAL_HELD = 16; // while the finger constraint is active
112+// The held sheet counts as settled when its non-rigid kinetic energy drops
113+// below settleKE, or when it has stopped falling: it has not halved over the
114+// last SETTLE_WINDOW seconds (or 3 / settleSigma0, whichever is longer).
115+export const SETTLE_WINDOW = 0.05; // [s]
116+export const SETTLE_RATIO = 0.5;
98117
99118 // ---------- Grid ----------
100119
@@ -250,6 +269,34 @@ export function dispersionFrequency(m: Material, g: Grid, err: number): number {
250269
251270 // ---------- Geometry of the excitation ----------
252271
272+/** One drumstick hit: time after the first impact [s] and stick speed [m/s]. */
273+export interface Hit {
274+ t: number;
275+ v: number;
276+}
277+
278+/**
279+ * Schedule of the drumstick hits. The first hit is always at t = 0. A series
280+ * has hitCount hits every hitInterval seconds, all at strikeVelocity for
281+ * "steady", rising from strikeVelocityMin to strikeVelocity for "crescendo",
282+ * and falling for "decrescendo".
283+ */
284+export function hitSchedule(p: PhysParams): Hit[] {
285+ if (p.hitSequence === "single") return [{ t: 0, v: p.strikeVelocity }];
286+ const count = Math.max(1, Math.round(p.hitCount));
287+ const dt = Math.max(0, p.hitInterval);
288+ const lo = Math.min(p.strikeVelocityMin, p.strikeVelocity), hi = p.strikeVelocity;
289+ const hits: Hit[] = [];
290+ for (let k = 0; k < count; k++) {
291+ const f = count > 1 ? k / (count - 1) : 1;
292+ let v = hi;
293+ if (p.hitSequence === "crescendo") v = lo + (hi - lo) * f;
294+ else if (p.hitSequence === "decrescendo") v = hi - (hi - lo) * f;
295+ hits.push({ t: k * dt, v });
296+ }
297+ return hits;
298+}
299+
253300 /** Strike point in sheet coordinates. */
254301 export function strikePoint(p: PhysParams): [number, number] {
255302 return [(p.strikeFracX - 0.5) * p.Lx, (p.strikeFracY - 0.5) * p.Ly];
src/selftest.tsmodified+3−1View file
@@ -23,6 +23,7 @@ const cases: { name: string; overrides: Partial<PhysParams> }[] = [
2323 { name: "rectangle, flat", overrides: { shape: "rectangle", bowShape: "none" } },
2424 { name: "rectangle, bowed across the width", overrides: { shape: "rectangle", bowShape: "width" } },
2525 { name: "ellipse, dish", overrides: { shape: "ellipse", bowShape: "dish" } },
26+ { name: "rectangle, bowed, in-plane mass x4", overrides: { shape: "rectangle", bowShape: "width", inPlaneMass: 4 } },
2627 ];
2728
2829 /** A smooth random field: a few sinusoids across the sheet plus a little white noise. */
@@ -75,9 +76,10 @@ export async function runSelfTest(device: GPUDevice): Promise<SelfTestCase[]> {
7576 const out = await sim.readState();
7677 sim.destroy();
7778
78- // Reference: a = elastic + losses, v' = v + dt a, w' = w + dt v'.
79+ // Reference: a = elastic (in-plane scaled) + losses, v' = v + dt a, w' = w + dt v'.
7980 const acc = fields(n);
8081 accelerations(g, m, x, acc, w0);
82+ for (let k = 0; k < n; k++) { acc.u[k] /= params.inPlaneMass; acc.v[k] /= params.inPlaneMass; }
8183 addDamping(g, params.sigma0, params.sigma1, params.sigmaIn, v, acc);
8284 const dt = sim.dt;
8385 const err = { u: 0, v: 0, w: 0 };
src/sim/cpuRef.tsmodified+5−1View file
@@ -203,8 +203,9 @@ export function potentialEnergy(g: Grid, m: Material, f: Fields, w0?: Float64Arr
203203 * Largest angular frequency of the linearized discrete system, by power
204204 * iteration on a small full grid with the same spacing. The extreme mode is a
205205 * grid-scale in-plane mode, so neither the sheet size nor its outline matters.
206+ * `inPlaneMass` is the factor by which the in-plane inertia is scaled.
206207 */
207-export function maxAngularFrequency(g: Grid, m: Material, nSmall = 12, iters = 300): number {
208+export function maxAngularFrequency(g: Grid, m: Material, inPlaneMass = 1, nSmall = 12, iters = 300): number {
208209 const gs: Grid = { nx: nSmall, ny: nSmall, hx: g.hx, hy: g.hy, n: nSmall * nSmall };
209210 const n = gs.n;
210211 const x = fields(n), y = fields(n);
@@ -226,6 +227,9 @@ export function maxAngularFrequency(g: Grid, m: Material, nSmall = 12, iters = 3
226227 x.u[k] *= eps / norm; x.v[k] *= eps / norm; x.w[k] *= eps / norm;
227228 }
228229 accelerations(gs, m, x, y);
230+ if (inPlaneMass !== 1) {
231+ for (let k = 0; k < n; k++) { y.u[k] /= inPlaneMass; y.v[k] /= inPlaneMass; }
232+ }
229233 let dot = 0;
230234 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];
231235 lambda = -dot / (eps * eps); // Rayleigh quotient of -M^{-1} K
src/sim/shaders/audio.wgslmodified+8−5View file
@@ -1,10 +1,13 @@
1-// Retarded-time sum over the plate for one output sample (Rayleigh integral).
1+// Retarded-time sum over the plate (Rayleigh integral), one workgroup per
2+// output sample of the batch, dispatched once at the end of the batch.
23 // lw[k] = (weight, delay in samples) for node k; delay is relative to the
34 // closest node so that the ring buffer only has to cover the spread.
45 // The ring buffer holds the transverse velocity at the sample instants; the
56 // acceleration is the first difference of the delayed velocity, which is the
67 // average acceleration over one sample interval (band-limited, and immune to
7-// the per-step rounding noise of the force evaluation).
8+// the per-step rounding noise of the force evaluation). The bound StepParams
9+// are those of the batch's first audio step; workgroup q handles the sample
10+// q slots later.
811
912 @group(0) @binding(0) var<uniform> P: Params;
1013 @group(0) @binding(1) var<uniform> SP: StepParams;
@@ -15,10 +18,10 @@
1518 var<workgroup> sh: array<f32, 256>;
1619
1720 @compute @workgroup_size(256)
18-fn main(@builtin(local_invocation_id) lid: vec3u) {
21+fn main(@builtin(local_invocation_id) lid: vec3u, @builtin(workgroup_id) wid: vec3u) {
1922 let t = lid.x;
2023 let H = i32(P.histLen);
21- let slot = SP.audioSlot;
24+ let slot = (SP.audioSlot + i32(wid.x)) % H;
2225 var acc = 0.0;
2326 for (var k = t; k < P.n; k += 256u) {
2427 let d = lw[k].y;
@@ -39,5 +42,5 @@ fn main(@builtin(local_invocation_id) lid: vec3u) {
3942 if (t < stride) { sh[t] += sh[t + stride]; }
4043 workgroupBarrier();
4144 }
42- if (t == 0u) { audioOut[u32(SP.audioOut)] = sh[0]; }
45+ if (t == 0u) { audioOut[u32(SP.audioOut) + wid.x] = sh[0]; }
4346 }
src/sim/shaders/common.wgslmodified+5−5View file
@@ -35,16 +35,16 @@ struct Params {
3535 patchJ1: u32,
3636 patchCX: f32, // mass-weighted center of the contact patch [m]
3737 patchCY: f32,
38- pad0: u32,
39- pad1: u32,
38+ inPlaneMass: f32, // factor by which the in-plane inertia is scaled (1: physical)
39+ invInPlaneMass: f32, // its reciprocal
4040 };
4141
4242 // Per-step parameters, addressed with a dynamic uniform offset.
4343 struct StepParams {
4444 step: u32,
45- flags: u32, // bit 0: push constraint active; bit 1: record its reaction; bit 2: apply the recorded reaction; bit 3: drumstick active
45+ flags: u32, // bit 0: push constraint active; bit 1: record its reaction; bit 2: apply the recorded reaction; bit 3: drumstick active; bit 4: re-arm the stick
4646 audioSlot: i32, // ring-buffer slot to record the velocity into (-1: none)
47- audioOut: i32, // output sample index within the batch
47+ audioOut: i32, // index of this sample within the batch's output
4848 pushW: f32, // prescribed patch displacement [m]
4949 pushV: f32, // prescribed patch velocity [m/s]
5050 sigma0: f32,
@@ -52,7 +52,7 @@ struct StepParams {
5252 releaseGain: f32, // fraction of the recorded reaction still applied while the finger lets go
5353 sigmaIn: f32, // additional frequency-independent loss of the in-plane motion
5454 parity: u32, // which of the two stick-state slots holds the current state
55- pad: u32,
55+ strikeV: f32, // stick speed of a re-armed hit [m/s] (bit 4)
5656 };
5757
5858 fn nodeX(i: i32) -> f32 { return (f32(i) - 0.5 * f32(P.nx - 1u)) * P.hx; }
src/sim/shaders/rigid.wgslmodified+6−3View file
@@ -1,6 +1,8 @@
11 // Rigid-body projection: (1) mass-weighted sums of the six rigid modes for
22 // displacement and velocity, kinetic energy and max |w| in one workgroup;
3-// (2) subtraction of the rigid components from every node.
3+// (2) subtraction of the rigid components from every node. The in-plane
4+// inertia may be scaled by P.inPlaneMass; the rigid coefficients are ratios
5+// and do not see the factor, the kinetic energy does.
46
57 @group(0) @binding(0) var<uniform> P: Params;
68 @group(0) @binding(1) var<storage, read_write> pos: array<f32>; // [u | v | w | wLo]
@@ -26,7 +28,7 @@ fn sums(@builtin(local_invocation_id) lid: vec3u) {
2628 s[3] += m * (x * p.y - y * p.x); s[4] += m * y * p.z; s[5] += m * (-x * p.z);
2729 s[6] += m * v.x; s[7] += m * v.y; s[8] += m * v.z;
2830 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);
31+ s[12] += 0.5 * m * (P.inPlaneMass * (v.x * v.x + v.y * v.y) + v.z * v.z);
3032 s[13] = max(s[13], abs(p.z));
3133 }
3234 for (var q = 0u; q < NS; q++) { sh[t][q] = s[q]; }
@@ -45,8 +47,9 @@ fn sums(@builtin(local_invocation_id) lid: vec3u) {
4547 stats[q + 6u] = sh[0][q + 6u] * inv[q];
4648 }
4749 // Kinetic energy of the non-rigid motion (the six rigid modes are mass-orthogonal).
50+ let ip = array<f32, 6>(P.inPlaneMass, P.inPlaneMass, 1.0, P.inPlaneMass, 1.0, 1.0);
4851 var keRigid = 0.0;
49- for (var q = 6u; q < 12u; q++) { keRigid += 0.5 * sh[0][q] * sh[0][q] * inv[q - 6u]; }
52+ for (var q = 6u; q < 12u; q++) { keRigid += 0.5 * ip[q - 6u] * sh[0][q] * sh[0][q] * inv[q - 6u]; }
5053 stats[12] = max(sh[0][12] - keRigid, 0.0);
5154 stats[13] = sh[0][13];
5255 }
src/sim/shaders/step.wgslmodified+94−49View file
@@ -3,11 +3,20 @@
33 // State layout (structure of arrays, n = nx * ny nodes):
44 // pos: [u | v | w | wLo] wLo carries the rounding error of the w update
55 // vel: [du | dv | dw]
6-// Each 16x16 workgroup loads its tile plus a 2-node halo of w and of the node
7-// weights into shared memory once, computes the membrane forces of the cells
8-// around the tile, and then evaluates the stencils from shared memory. Nodes
9-// outside the sheet (weight 0) are skipped: a second difference exists only
10-// where the node and both neighbors belong to the sheet, and a cell
6+// A workgroup of 16x16 threads advances a 12x12 tile of nodes in three
7+// phases separated by barriers. (1) Every thread loads one entry of the tile
8+// plus its 2-node halo (16x16 values of w, u, v, w0, w_t and the node
9+// weights) into shared memory. (2) The threads compute, once each and again
10+// into shared memory, the membrane forces of the 13x13 cells around the tile
11+// and the weighted bending moments and mixed differences that the bending
12+// stencil gathers. (3) Each of the 144 tile nodes evaluates its update from
13+// shared memory alone. On grids of a few thousand nodes the kernel is bound
14+// by the latency of one workgroup rather than by throughput, so the shared
15+// phases, which replace a tenfold per-node recomputation of the stencil
16+// ingredients, are what make a step cheap.
17+//
18+// Nodes outside the sheet (weight 0) are skipped: a second difference exists
19+// only where the node and both neighbors belong to the sheet, and a cell
1120 // contributes only when its four corners do. That reproduces the free edge of
1221 // the rectangle and gives a staircase free edge for an ellipse.
1322
@@ -23,15 +32,24 @@
2332 @group(0) @binding(7) var<storage, read_write> react: array<f32>;
2433 @group(0) @binding(8) var<storage, read> w0: array<f32>; // stress-free initial shape
2534
26-const TILE = 16;
35+const TILE = 12; // nodes per workgroup along each axis (must match simulation.ts)
2736 const HALO = 2;
28-const TW = 20; // TILE + 2 * HALO
29-var<workgroup> tw: array<f32, 400>;
30-// Node weights of the same tile in quarters (0..4), packed four to a word.
31-var<workgroup> tn: array<u32, 100>;
32-// Membrane force contributions of the 17x17 cells around the tile, 4 nodes x 3 components each.
33-const CT = 17;
34-var<workgroup> cellF: array<f32, 3468>;
37+const TW = 16; // TILE + 2 * HALO: the loaded tile, one thread per entry
38+const CT = 13; // TILE + 1: cells adjacent to the tile's nodes
39+const MT = 14; // TILE + 2: nodes at which the bending moments are needed
40+// Tile with halo, local index (j - oy) * TW + (i - ox).
41+var<workgroup> tw: array<f32, 256>; // w
42+var<workgroup> tu: array<f32, 256>; // u
43+var<workgroup> tv: array<f32, 256>; // v
44+var<workgroup> tw0: array<f32, 256>; // stress-free shape
45+var<workgroup> tdw: array<f32, 256>; // w_t
46+var<workgroup> tn: array<u32, 64>; // node weights in quarters (0..4), four to a word
47+// Membrane force contributions of the CT x CT cells, 4 nodes x 3 components each.
48+var<workgroup> cellF: array<f32, 2028>;
49+// Weighted bending moments at the MT x MT inner nodes and mixed differences at the CT x CT cells.
50+var<workgroup> ma: array<f32, 196>;
51+var<workgroup> mb: array<f32, 196>;
52+var<workgroup> cc: array<f32, 169>;
3553
3654 // Tile origin in global node coordinates (set once per invocation).
3755 var<private> ox: i32;
@@ -42,16 +60,18 @@ var<private> ihxy: f32;
4260
4361 fn id(i: i32, j: i32) -> u32 { return u32(j) * P.nx + u32(i); }
4462 fn lidx(i: i32, j: i32) -> u32 { return u32(j - oy) * u32(TW) + u32(i - ox); }
63+fn midx(i: i32, j: i32) -> u32 { return u32(j - oy - 1) * u32(MT) + u32(i - ox - 1); }
64+fn cidx(ci: i32, cj: i32) -> u32 { return u32(cj - oy - 1) * u32(CT) + u32(ci - ox - 1); }
4565 // Tile lookups with global indices; the tile covers [ox, ox + TW) x [oy, oy + TW).
4666 fn W(i: i32, j: i32) -> f32 { return tw[lidx(i, j)]; }
67+fn U(i: i32, j: i32) -> f32 { return tu[lidx(i, j)]; }
68+fn V(i: i32, j: i32) -> f32 { return tv[lidx(i, j)]; }
69+fn W0(i: i32, j: i32) -> f32 { return tw0[lidx(i, j)]; }
70+fn DW(i: i32, j: i32) -> f32 { return tdw[lidx(i, j)]; }
4771 fn NWq(i: i32, j: i32) -> u32 { let l = lidx(i, j); return (tn[l >> 2u] >> ((l & 3u) * 8u)) & 0xFFu; }
4872 fn NW(i: i32, j: i32) -> f32 { return 0.25 * f32(NWq(i, j)); }
4973 fn INS(i: i32, j: i32) -> bool { return NWq(i, j) != 0u; }
5074 fn insCell(ci: i32, cj: i32) -> bool { return INS(ci, cj) && INS(ci + 1, cj) && INS(ci + 1, cj + 1) && INS(ci, cj + 1); }
51-fn U(i: i32, j: i32) -> f32 { return posIn[id(i, j)]; }
52-fn V(i: i32, j: i32) -> f32 { return posIn[P.n + id(i, j)]; }
53-fn DW(i: i32, j: i32) -> f32 { return velIn[2u * P.n + id(i, j)]; }
54-fn W0(i: i32, j: i32) -> f32 { return w0[id(i, j)]; }
5575
5676 // Graph Laplacian of the transverse velocity at node (i, j), over neighbors in the sheet.
5777 fn lapDW(i: i32, j: i32) -> f32 {
@@ -137,45 +157,58 @@ fn inContactPatch(i: i32, j: i32) -> bool {
137157
138158 // Stored membrane force on local node `a` of cell (ci, cj).
139159 fn cellForce(ci: i32, cj: i32, a: u32) -> vec3f {
140- let lc = u32(cj - oy - 1) * u32(CT) + u32(ci - ox - 1);
141- let base = lc * 12u + a * 3u;
160+ let base = cidx(ci, cj) * 12u + a * 3u;
142161 return vec3f(cellF[base], cellF[base + 1u], cellF[base + 2u]);
143162 }
144163
145164 @compute @workgroup_size(16, 16)
146-fn main(@builtin(global_invocation_id) gid: vec3u, @builtin(local_invocation_id) lid: vec3u, @builtin(workgroup_id) wid: vec3u) {
165+fn main(@builtin(local_invocation_id) lid: vec3u, @builtin(workgroup_id) wid: vec3u) {
147166 let n = P.n;
148167 ox = i32(wid.x) * TILE - HALO;
149168 oy = i32(wid.y) * TILE - HALO;
150169 ihx2 = 1.0 / (P.hx * P.hx); ihy2 = 1.0 / (P.hy * P.hy); ihxy = 1.0 / (P.hx * P.hy);
151- // Cooperative load of the w and weight tiles (400 entries each, 256 threads).
152170 let lt = lid.y * 16u + lid.x;
153- for (var e = lt; e < 400u; e += 256u) {
154- let gi = ox + i32(e % u32(TW)); let gj = oy + i32(e / u32(TW));
155- var ww = 0.0;
156- if (gi >= 0 && gi < i32(P.nx) && gj >= 0 && gj < i32(P.ny)) { ww = posIn[2u * n + id(gi, gj)]; }
157- tw[e] = ww;
158- }
159- if (lt < 100u) {
160- var word = 0u;
161- for (var q = 0u; q < 4u; q++) {
162- let e = lt * 4u + q;
163- let gi = ox + i32(e % u32(TW)); let gj = oy + i32(e / u32(TW));
164- var wq = 0u;
165- if (gi >= 0 && gi < i32(P.nx) && gj >= 0 && gj < i32(P.ny)) { wq = u32(round(4.0 * nw[id(gi, gj)])); }
166- word |= wq << (q * 8u);
171+
172+ // Phase 1: one entry of the tile with halo per thread.
173+ {
174+ let gi = ox + i32(lid.x); let gj = oy + i32(lid.y);
175+ var ww = 0.0; var uu = 0.0; var vv = 0.0; var z0 = 0.0; var dw = 0.0;
176+ if (gi >= 0 && gi < i32(P.nx) && gj >= 0 && gj < i32(P.ny)) {
177+ let k = id(gi, gj);
178+ ww = posIn[2u * n + k]; uu = posIn[k]; vv = posIn[n + k]; z0 = w0[k]; dw = velIn[2u * n + k];
179+ }
180+ tw[lt] = ww; tu[lt] = uu; tv[lt] = vv; tw0[lt] = z0; tdw[lt] = dw;
181+ if (lt < 64u) {
182+ var word = 0u;
183+ for (var q = 0u; q < 4u; q++) {
184+ let e = lt * 4u + q;
185+ let ei = ox + i32(e % u32(TW)); let ej = oy + i32(e / u32(TW));
186+ var wq = 0u;
187+ if (ei >= 0 && ei < i32(P.nx) && ej >= 0 && ej < i32(P.ny)) { wq = u32(round(4.0 * nw[id(ei, ej)])); }
188+ word |= wq << (q * 8u);
189+ }
190+ tn[lt] = word;
167191 }
168- tn[lt] = word;
169192 }
170193 workgroupBarrier();
171- // Membrane forces of the 17x17 cells adjacent to the tile's nodes.
172- for (var e = lt; e < 289u; e += 256u) {
173- let lci = i32(e % u32(CT)); let lcj = i32(e / u32(CT));
174- cellForces(ox + 1 + lci, oy + 1 + lcj, e);
194+
195+ // Phase 2: membrane forces of the CT x CT cells, bending moments at the MT x MT
196+ // inner nodes, mixed differences at the cells; one item per thread each.
197+ if (lt < 169u) {
198+ let ci = ox + 1 + i32(lt % u32(CT)); let cj = oy + 1 + i32(lt / u32(CT));
199+ cellForces(ci, cj, lt);
200+ cc[lt] = Cc(ci, cj);
201+ }
202+ if (lt < 196u) {
203+ let mi = ox + 1 + i32(lt % u32(MT)); let mj = oy + 1 + i32(lt / u32(MT));
204+ ma[lt] = Ma(mi, mj);
205+ mb[lt] = Mb(mi, mj);
175206 }
176207 workgroupBarrier();
177208
178- let i = i32(gid.x); let j = i32(gid.y);
209+ // Phase 3: the update of the TILE x TILE nodes (no barriers past this point).
210+ if (lt >= 144u) { return; }
211+ let i = ox + HALO + i32(lt % u32(TILE)); let j = oy + HALO + i32(lt / u32(TILE));
179212 if (i >= i32(P.nx) || j >= i32(P.ny)) { return; }
180213 let k = id(i, j);
181214 let weight = NW(i, j);
@@ -188,9 +221,9 @@ fn main(@builtin(global_invocation_id) gid: vec3u, @builtin(local_invocation_id)
188221 }
189222
190223 // Elastic forces: bending from the discrete Kirchhoff energy, membrane from the cells.
191- let bend = (Ma(i - 1, j) - 2.0 * Ma(i, j) + Ma(i + 1, j)) * ihx2
192- + (Mb(i, j - 1) - 2.0 * Mb(i, j) + Mb(i, j + 1)) * ihy2
193- + 2.0 * (1.0 - P.nu) * (Cc(i, j) - Cc(i - 1, j) - Cc(i, j - 1) + Cc(i - 1, j - 1)) * ihxy;
224+ let bend = (ma[midx(i - 1, j)] - 2.0 * ma[midx(i, j)] + ma[midx(i + 1, j)]) * ihx2
225+ + (mb[midx(i, j - 1)] - 2.0 * mb[midx(i, j)] + mb[midx(i, j + 1)]) * ihy2
226+ + 2.0 * (1.0 - P.nu) * (cc[cidx(i, j)] - cc[cidx(i - 1, j)] - cc[cidx(i, j - 1)] + cc[cidx(i - 1, j - 1)]) * ihxy;
194227 var F = vec3f(0.0, 0.0, -P.D * P.hx * P.hy * bend);
195228 F += cellForce(i, j, 0u);
196229 F += cellForce(i - 1, j, 1u);
@@ -198,6 +231,8 @@ fn main(@builtin(global_invocation_id) gid: vec3u, @builtin(local_invocation_id)
198231 F += cellForce(i, j - 1, 3u);
199232 let mass = P.rhoh * P.hx * P.hy * weight;
200233 var acc = F * (1.0 / mass);
234+ // Optional scaling of the in-plane inertia (lengthens the stable time step).
235+ acc.x *= P.invInPlaneMass; acc.y *= P.invInPlaneMass;
201236
202237 // Losses: frequency-independent on all components, the graph Laplacian of the
203238 // transverse velocity (~k^2), and an extra frequency-independent in-plane loss.
@@ -238,11 +273,21 @@ fn main(@builtin(global_invocation_id) gid: vec3u, @builtin(local_invocation_id)
238273 acc.z -= Fc / P.patchMass;
239274 if (k == P.strikeNode) {
240275 let sOut = n + 4u * (1u - SP.parity);
241- let svn = sv + P.dt * Fc / P.stickMass;
242- react[sOut] = sx + P.dt * svn;
243- react[sOut + 1u] = svn;
244- react[sOut + 2u] = Fc;
245- react[sOut + 3u] = max(sMax, Fc);
276+ if ((SP.flags & 16u) != 0u) {
277+ // Next hit of a series: a fresh stick arrives at the patch, wherever
278+ // it is at this instant, with the scheduled speed; the previous
279+ // stick's rebound is discarded.
280+ react[sOut] = wc;
281+ react[sOut + 1u] = -SP.strikeV;
282+ react[sOut + 2u] = 0.0;
283+ react[sOut + 3u] = sMax;
284+ } else {
285+ let svn = sv + P.dt * Fc / P.stickMass;
286+ react[sOut] = sx + P.dt * svn;
287+ react[sOut + 1u] = svn;
288+ react[sOut + 2u] = Fc;
289+ react[sOut + 3u] = max(sMax, Fc);
290+ }
246291 }
247292 }
248293 }
src/sim/simulation.tsmodified+149−48View file
@@ -3,9 +3,12 @@
33 // The simulation advances in batches of one movie frame (1/fps seconds). A
44 // batch encodes one compute dispatch per time step (with a per-step uniform
55 // addressed by dynamic offset), a rigid-body projection every RB_INTERVAL
6-// steps, one retarded-time listener sum per audio sample, and one movie-frame
7-// extraction at the end. Audio, frame, and statistics are read back after
8-// every batch.
6+// steps, one retarded-time listener sum over the batch's audio samples, and
7+// one movie-frame extraction at the end. Audio, frame, and statistics are
8+// read back after every batch. One batch is kept encoded ahead of the one
9+// being read back, so the GPU does not wait for the readback or for the
10+// encoding; a decision that depends on a readback (the release of the held
11+// sheet) therefore takes effect one batch later.
912
1013 import {
1114 type PhysParams,
@@ -17,14 +20,19 @@ import {
1720 strikePoint,
1821 pushPoint,
1922 contactPatch,
23+ hitSchedule,
24+ type Hit,
2025 nodeX,
2126 nodeY,
2227 nodeWeight,
2328 nodeWeights,
24- HIST_LEN,
29+ CFL,
30+ HIST_MARGIN,
2531 RB_INTERVAL,
2632 RB_INTERVAL_HELD,
2733 RHO_AIR,
34+ SETTLE_RATIO,
35+ SETTLE_WINDOW,
2836 SPEED_OF_SOUND,
2937 } from "../params.ts";
3038 import { maxAngularFrequency } from "./cpuRef.ts";
@@ -56,12 +64,28 @@ export interface SimStatus {
5664 delayClamped: boolean;
5765 /** Drumstick state (strike excitation only). */
5866 stick: StickState | null;
67+ /** Hits delivered so far, of the scheduled total (strike excitation only). */
68+ hits: { done: number; total: number } | null;
5969 }
6070
6171 const PARAMS_LEN = 36; // floats in the Params uniform
6272 const STEP_LEN = 12; // floats in the StepParams uniform
6373 const STATS_LEN = 16;
6474 const STICK_LEN = 8;
75+/** Nodes per workgroup along each axis in step.wgsl. */
76+const TILE = 12;
77+
78+/** A batch that has been encoded and submitted but not read back. */
79+interface Pending {
80+ staging: GPUBuffer;
81+ steps: number;
82+ nAudio: number;
83+ recording: boolean;
84+ /** Preparation clock at the end of the batch [s]. */
85+ prepTimeEnd: number;
86+ /** Which of the two stick-state slots the batch's last step wrote. */
87+ stickParity: number;
88+}
6589
6690 export class Simulation {
6791 readonly grid: Grid;
@@ -72,6 +96,8 @@ export class Simulation {
7296 readonly samplesPerBatch: number;
7397 readonly stepsPerBatch: number;
7498 readonly omegaMax: number;
99+ /** Ring buffer length of the listener sum [samples]. */
100+ readonly histLen: number;
75101
76102 /** Listener pressure [Pa] at fs, from the release or the impact. */
77103 audio: Float32Array = new Float32Array(0);
@@ -86,12 +112,22 @@ export class Simulation {
86112 peakPressure = 0;
87113 delayClamped = false;
88114
89- /** Steps since the start of the current clock (preparation or recording). */
115+ /** Steps enqueued since the start of the current clock (preparation or recording). */
90116 private stepCount = 0;
117+ /** Preparation time enqueued [s]. */
91118 private prepTime = 0;
92119 private released = false;
93120 private lastPrepFrame: Float32Array<ArrayBuffer> | null = null;
94121 private stick: StickState | null = null;
122+ /** Drumstick hit schedule (strike excitation only) and the index of the next hit to arm. */
123+ readonly hits: Hit[];
124+ private hitSteps: number[];
125+ private nextHit = 1;
126+ private inFlight: Pending | null = null;
127+ /** The held sheet has settled; release once the batch in flight has been read back. */
128+ private releaseRequested = false;
129+ /** Kinetic energy after each batch of the hold, for the plateau criterion. */
130+ private keHistory: { t: number; ke: number }[] = [];
95131
96132 private device: GPUDevice;
97133 private params: PhysParams;
@@ -141,12 +177,13 @@ export class Simulation {
141177 this.mat = m;
142178
143179 // Stable time step, rounded so that an integer number of steps make one audio sample.
144- this.omegaMax = maxAngularFrequency(g, m);
145- const dtMax = (0.85 * 2) / this.omegaMax;
180+ this.omegaMax = maxAngularFrequency(g, m, params.inPlaneMass);
181+ const dtMax = (CFL * 2) / this.omegaMax;
146182 this.stepsPerSample = Math.ceil(1 / (params.fs * dtMax));
147183 this.dt = 1 / (params.fs * this.stepsPerSample);
148184 this.samplesPerBatch = Math.max(1, Math.round(params.fs / params.fps));
149185 this.stepsPerBatch = this.samplesPerBatch * this.stepsPerSample;
186+ this.histLen = this.samplesPerBatch + HIST_MARGIN;
150187
151188 const n = g.n;
152189 const f32 = (size: number, usage = GPUBufferUsage.STORAGE | GPUBufferUsage.COPY_SRC | GPUBufferUsage.COPY_DST) =>
@@ -154,7 +191,7 @@ export class Simulation {
154191
155192 this.pos = [f32(4 * n), f32(4 * n)]; // [u | v | w | wLo]
156193 this.vel = [f32(3 * n), f32(3 * n)]; // [du | dv | dw]
157- this.hist = f32(HIST_LEN * n);
194+ this.hist = f32(this.histLen * n);
158195 this.react = f32(n + STICK_LEN);
159196 this.w0 = f32(n);
160197 device.queue.writeBuffer(this.w0, 0, initialShape(params, g));
@@ -180,17 +217,21 @@ export class Simulation {
180217 Iy += mk * nodeX(g, i) ** 2;
181218 }
182219 }
183- // Drumstick contact patch and initial stick state (position 0, moving toward the sheet).
220+ // Drumstick contact patch, hit schedule (in steps, strictly increasing),
221+ // and initial stick state (position 0, moving toward the sheet at the speed of the first hit).
184222 const cp = contactPatch(params, g);
223+ this.hits = hitSchedule(params);
224+ this.hitSteps = [];
225+ for (const h of this.hits) this.hitSteps.push(Math.max(this.hitSteps.length === 0 ? 0 : this.hitSteps[this.hitSteps.length - 1] + 1, Math.round(h.t / this.dt)));
185226 const stick = new Float32Array(STICK_LEN);
186- stick[1] = -params.strikeVelocity; stick[5] = -params.strikeVelocity;
227+ stick[1] = -this.hits[0].v; stick[5] = -this.hits[0].v;
187228 device.queue.writeBuffer(this.react, n * 4, stick);
188229
189230 const [px, py] = pushPoint(params);
190231 const [sx, sy] = strikePoint(params);
191232 const pf = new Float32Array(PARAMS_LEN);
192233 const pu = new Uint32Array(pf.buffer);
193- pu[0] = g.nx; pu[1] = g.ny; pu[2] = n; pu[3] = HIST_LEN;
234+ pu[0] = g.nx; pu[1] = g.ny; pu[2] = n; pu[3] = this.histLen;
194235 pf[4] = g.hx; pf[5] = g.hy; pf[6] = this.dt; pf[7] = m.nu;
195236 pf[8] = m.D; pf[9] = m.C; pf[10] = m.G; pf[11] = m.rhoh;
196237 pf[12] = px; pf[13] = py; pf[14] = params.fingerRadius ** 2; pf[15] = mass;
@@ -199,6 +240,7 @@ export class Simulation {
199240 pu[24] = cp.nearest; pf[25] = params.contactStiffness; pf[26] = params.contactExponent; pf[27] = params.stickMass;
200241 pu[28] = cp.box[0]; pu[29] = cp.box[1]; pu[30] = cp.box[2]; pu[31] = cp.box[3];
201242 pf[32] = cp.center[0]; pf[33] = cp.center[1];
243+ pf[34] = params.inPlaneMass; pf[35] = 1 / params.inPlaneMass;
202244 this.paramsBuf = device.createBuffer({ size: pf.byteLength, usage: GPUBufferUsage.UNIFORM | GPUBufferUsage.COPY_DST });
203245 device.queue.writeBuffer(this.paramsBuf, 0, pf);
204246
@@ -347,8 +389,8 @@ export class Simulation {
347389 this.delayClamped = false;
348390 for (let k = 0; k < g.n; k++) {
349391 let d = data[2 * k + 1] - dmin;
350- if (d > HIST_LEN - 3) {
351- d = HIST_LEN - 3;
392+ if (d > HIST_MARGIN - 3) {
393+ d = HIST_MARGIN - 3;
352394 this.delayClamped = true;
353395 }
354396 data[2 * k + 1] = d;
@@ -385,7 +427,9 @@ export class Simulation {
385427 }
386428
387429 status(): SimStatus {
388- const t = this.phase === "push" || this.phase === "hold" ? this.prepTime - this.prepDuration() : this.stepCount * this.dt;
430+ // Report the time of the state that has been read back, not the one enqueued.
431+ const ahead = this.inFlight ? this.inFlight.steps * this.dt : 0;
432+ const t = this.phase === "push" || this.phase === "hold" ? this.prepTime - ahead - this.prepDuration() : this.stepCount * this.dt - ahead;
389433 return {
390434 phase: this.phase,
391435 t,
@@ -397,6 +441,7 @@ export class Simulation {
397441 peakPressure: this.peakPressure,
398442 delayClamped: this.delayClamped,
399443 stick: this.stick,
444+ hits: this.params.excitation === "strike" ? { done: this.hitSteps.filter((h) => h * this.dt <= t).length, total: this.hits.length } : null,
400445 };
401446 }
402447
@@ -413,6 +458,10 @@ export class Simulation {
413458 pass.end();
414459 }
415460
461+ private stepWorkgroups(): [number, number] {
462+ return [Math.ceil(this.grid.nx / TILE), Math.ceil(this.grid.ny / TILE)];
463+ }
464+
416465 /** The recorded listener pressure so far [Pa]. */
417466 recordedAudio(): Float32Array {
418467 return this.audio.subarray(0, this.samples);
@@ -431,7 +480,7 @@ export class Simulation {
431480 }
432481
433482 /** Writes the per-step uniform for step s of a batch. */
434- private writeStepParams(s: number, fields: { flags: number; slot: number; outIdx: number; pushW: number; pushV: number; sigma0: number; releaseGain: number }): void {
483+ private writeStepParams(s: number, fields: { flags: number; slot: number; outIdx: number; pushW: number; pushV: number; sigma0: number; releaseGain: number; strikeV?: number }): void {
435484 const o = (s * this.stepStride) / 4;
436485 const u32 = new Uint32Array(this.stepData);
437486 const i32 = new Int32Array(this.stepData);
@@ -447,27 +496,47 @@ export class Simulation {
447496 f32[o + 8] = fields.releaseGain;
448497 f32[o + 9] = this.params.sigmaIn;
449498 u32[o + 10] = (this.totalSteps + s) & 1;
450- u32[o + 11] = 0;
499+ f32[o + 11] = fields.strikeV ?? 0;
451500 }
452501
453- /** Runs one batch (one movie frame of time). Resolves when its results have been read back. */
502+ /**
503+ * Runs one batch (one movie frame of time). Resolves when its results have
504+ * been read back. The following batch is encoded before this one is awaited,
505+ * unless its kind depends on this one's result.
506+ */
454507 async runBatch(): Promise<void> {
455508 if (this.phase === "done") return;
456- const p = this.params;
457- const recording = this.phase === "recording";
458- const stepsLeft = recording ? (this.targetSamples() - this.samples) * this.stepsPerSample : this.stepsPerBatch;
459- const steps = Math.min(this.stepsPerBatch, Math.max(0, stepsLeft));
460- if (steps === 0) {
509+ const cur = this.inFlight ?? this.encodeBatch();
510+ this.inFlight = null;
511+ if (!cur) {
461512 this.phase = "done";
462513 return;
463514 }
515+ if (this.canEncodeAhead()) this.inFlight = this.encodeBatch();
516+ await this.collect(cur);
517+ }
518+
519+ /** True when the batch after the one in flight is certainly needed and of a known kind. */
520+ private canEncodeAhead(): boolean {
521+ if (this.phase === "recording") return this.targetSamples() * this.stepsPerSample - this.stepCount > 0;
522+ return !this.releaseRequested && this.prepTime < this.prepDuration();
523+ }
524+
525+ /** Encodes and submits one batch; null when no steps remain. */
526+ private encodeBatch(): Pending | null {
527+ const p = this.params;
528+ const recording = this.phase === "recording";
529+ const stepsLeft = recording ? this.targetSamples() * this.stepsPerSample - this.stepCount : this.stepsPerBatch;
530+ const steps = Math.min(this.stepsPerBatch, Math.max(0, stepsLeft));
531+ if (steps === 0) return null;
464532
465533 // Per-step parameters: the push ramp and settle damping before the release,
466- // the release ramp or the drumstick and the audio bookkeeping afterwards.
534+ // the release ramp or the drumstick (re-armed at each scheduled hit) and
535+ // the audio bookkeeping afterwards.
467536 const audioSteps: number[] = [];
468537 for (let s = 0; s < steps; s++) {
469538 const local = this.stepCount + s;
470- let flags = 0, pushW = 0, pushV = 0, sigma0 = p.sigma0, slot = -1, outIdx = -1, releaseGain = 0;
539+ let flags = 0, pushW = 0, pushV = 0, sigma0 = p.sigma0, slot = -1, outIdx = -1, releaseGain = 0, strikeV = 0;
471540 if (!recording) {
472541 const t = this.prepTime + (s + 1) * this.dt; // state after this step
473542 const tau = Math.min(1, t / p.rampTime);
@@ -481,18 +550,23 @@ export class Simulation {
481550 const t = local * this.dt;
482551 if (p.excitation === "strike") {
483552 flags = 8;
553+ if (this.nextHit < this.hitSteps.length && this.hitSteps[this.nextHit] === local) {
554+ flags |= 16;
555+ strikeV = this.hits[this.nextHit].v;
556+ this.nextHit++;
557+ }
484558 } else if (p.releaseTime > 0 && t < p.releaseTime) {
485559 const tau = t / p.releaseTime;
486560 releaseGain = 1 - tau * tau * (3 - 2 * tau);
487561 flags = 4;
488562 }
489563 if (local % this.stepsPerSample === 0) {
490- slot = (local / this.stepsPerSample) % HIST_LEN;
564+ slot = (local / this.stepsPerSample) % this.histLen;
491565 outIdx = audioSteps.length;
492566 audioSteps.push(s);
493567 }
494568 }
495- this.writeStepParams(s, { flags, slot, outIdx, pushW, pushV, sigma0, releaseGain });
569+ this.writeStepParams(s, { flags, slot, outIdx, pushW, pushV, sigma0, releaseGain, strikeV });
496570 }
497571 this.device.queue.writeBuffer(this.stepBuf, 0, this.stepData, 0, steps * this.stepStride);
498572
@@ -500,10 +574,9 @@ export class Simulation {
500574 const CHUNK = 1024;
501575 let enc = this.device.createCommandEncoder();
502576 let pass = enc.beginComputePass();
503- const wgX = Math.ceil(this.grid.nx / 16), wgY = Math.ceil(this.grid.ny / 16);
577+ const [wgX, wgY] = this.stepWorkgroups();
504578 const wgN = Math.ceil(this.grid.n / 256);
505579 let cur = this.cur;
506- let audioPtr = 0;
507580 // While held, the projection defines the inertia-relief holding and runs
508581 // every RB_INTERVAL_HELD steps. While ringing it is optional; the sums
509582 // still run once per batch for the statistics.
@@ -537,13 +610,14 @@ export class Simulation {
537610 }
538611 stepBound = false;
539612 }
540- if (audioPtr < audioSteps.length && audioSteps[audioPtr] === s) {
541- pass.setPipeline(this.audioPipe);
542- pass.setBindGroup(0, this.audioBG, [s * this.stepStride]);
543- pass.dispatchWorkgroups(1);
544- audioPtr++;
545- stepBound = false;
546- }
613+ }
614+ // The listener sum for every audio sample of the batch, one workgroup each;
615+ // the ring buffer is long enough that no slot of the batch has been reused.
616+ const nAudio = audioSteps.length;
617+ if (nAudio > 0) {
618+ pass.setPipeline(this.audioPipe);
619+ pass.setBindGroup(0, this.audioBG, [audioSteps[0] * this.stepStride]);
620+ pass.dispatchWorkgroups(nAudio);
547621 }
548622 pass.setPipeline(this.extractPipe);
549623 pass.setBindGroup(0, this.extractBG[cur]);
@@ -551,7 +625,6 @@ export class Simulation {
551625 pass.end();
552626 const staging = this.staging[this.stagingIdx];
553627 this.stagingIdx = 1 - this.stagingIdx;
554- const nAudio = audioSteps.length;
555628 const n = this.grid.n;
556629 enc.copyBufferToBuffer(this.audioOut, 0, staging, 0, Math.max(4, nAudio * 4));
557630 enc.copyBufferToBuffer(this.frameBuf, 0, staging, this.samplesPerBatch * 4, n * 4);
@@ -562,38 +635,65 @@ export class Simulation {
562635 this.totalSteps += steps;
563636 this.stepCount += steps;
564637 if (!recording) this.prepTime += steps * this.dt;
638+ return { staging, steps, nAudio, recording, prepTimeEnd: this.prepTime, stickParity: this.totalSteps & 1 };
639+ }
565640
566- await staging.mapAsync(GPUMapMode.READ);
567- const view = new Float32Array(staging.getMappedRange());
641+ /** Reads back a batch and applies its results and the phase transitions they imply. */
642+ private async collect(b: Pending): Promise<void> {
643+ const p = this.params;
644+ const n = this.grid.n;
645+ await b.staging.mapAsync(GPUMapMode.READ);
646+ const view = new Float32Array(b.staging.getMappedRange());
568647 const stats = view.subarray(this.samplesPerBatch + n, this.samplesPerBatch + n + STATS_LEN);
569648 this.ke = stats[12];
570649 this.wmax = stats[13];
571650 if (p.excitation === "strike") {
572651 // The slot written by the last step is the one the next step reads.
573- const base = this.samplesPerBatch + n + STATS_LEN + 4 * (this.totalSteps & 1);
652+ const base = this.samplesPerBatch + n + STATS_LEN + 4 * b.stickParity;
574653 this.stick = { x: view[base], v: view[base + 1], force: view[base + 2], maxForce: view[base + 3] };
575654 }
576655 const frame = view.slice(this.samplesPerBatch, this.samplesPerBatch + n);
577- if (recording) {
578- this.appendAudio(view.slice(0, nAudio));
656+ if (b.recording) {
657+ this.appendAudio(view.slice(0, b.nAudio));
579658 this.frames.push(frame);
580659 } else {
581660 this.lastPrepFrame = frame;
582661 }
583- staging.unmap();
584-
585- if (!recording) {
586- if (this.phase === "push" && this.prepTime >= p.rampTime) this.phase = "hold";
587- const settled = this.phase === "hold" && this.ke < p.settleKE;
588- if (settled || this.prepTime >= this.prepDuration()) this.release();
662+ b.staging.unmap();
663+
664+ if (!b.recording) {
665+ if (this.phase === "push" && b.prepTimeEnd >= p.rampTime) this.phase = "hold";
666+ if (this.phase === "hold" && !this.releaseRequested && this.settled(b.prepTimeEnd)) this.releaseRequested = true;
667+ if (b.prepTimeEnd >= this.prepDuration()) this.releaseRequested = true;
668+ // Release only once every preparation batch has been read back, so that
669+ // frame 0 is the last held state.
670+ if (this.releaseRequested && !this.inFlight) this.release();
589671 } else if (this.samples >= this.targetSamples()) {
590672 this.phase = "done";
591673 }
592674 }
593675
676+ /**
677+ * Settle criterion of the held sheet at hold time t: the non-rigid kinetic
678+ * energy is below settleKE, or it has stopped falling (it has not halved
679+ * over the last window). The second case catches the rounding floor, which
680+ * can lie above any fixed threshold.
681+ */
682+ private settled(t: number): boolean {
683+ const p = this.params;
684+ this.keHistory.push({ t, ke: this.ke });
685+ if (this.ke < p.settleKE) return true;
686+ const window = Math.max(SETTLE_WINDOW, 3 / p.settleSigma0);
687+ let ref: { t: number; ke: number } | null = null;
688+ for (const h of this.keHistory) if (h.t <= t - window) ref = h;
689+ return ref !== null && this.ke >= SETTLE_RATIO * ref.ke;
690+ }
691+
594692 private release(): void {
595693 this.phase = "recording";
596694 this.released = true;
695+ this.releaseRequested = false;
696+ this.keHistory = [];
597697 this.stepCount = 0;
598698 this.totalSteps = 0;
599699 // The finger force and the distributed holding reaction relax together:
@@ -622,7 +722,7 @@ export class Simulation {
622722 const pass = enc.beginComputePass();
623723 pass.setPipeline(this.stepPipe);
624724 pass.setBindGroup(0, this.stepBG[this.cur], [0]);
625- pass.dispatchWorkgroups(Math.ceil(this.grid.nx / 16), Math.ceil(this.grid.ny / 16));
725+ pass.dispatchWorkgroups(...this.stepWorkgroups());
626726 pass.end();
627727 this.device.queue.submit([enc.finish()]);
628728 this.cur = 1 - this.cur;
@@ -654,6 +754,7 @@ export class Simulation {
654754 }
655755
656756 destroy(): void {
757+ this.inFlight = null;
657758 for (const b of [...this.pos, ...this.vel, this.hist, this.react, this.w0, this.nw, this.lw, this.audioOut, this.frameBuf, this.displayW, this.statsBuf, this.paramsBuf, this.stepBuf, ...this.staging]) b.destroy();
658759 }
659760
src/style.cssmodified+2−0View file
@@ -58,3 +58,5 @@ button.primary:hover { background: #f8c453; }
5858 .seg button.on { background: var(--accent2); border-color: var(--accent2); color: white; }
5959 .field[hidden] { display: none; }
6060 .field select.unit-select { padding: 3px 2px; font-size: 11px; }
61+.cacherow { display: flex; justify-content: space-between; align-items: center; gap: 6px; margin: 3px 0; font-size: 12px; color: var(--muted); }
62+.cacherow span { overflow-wrap: anywhere; }
src/ui/form.tsmodified+85−24View file
@@ -5,6 +5,7 @@ import {
55 type BowShape,
66 type Excitation,
77 type Grid,
8+ hitSchedule,
89 type PhysParams,
910 type PlateShape,
1011 countInside,
@@ -12,7 +13,9 @@ import {
1213 makeGrid,
1314 makeMaterial,
1415 t60,
16+ CFL,
1517 } from "../params.ts";
18+import { maxAngularFrequency } from "../sim/cpuRef.ts";
1619
1720 export type SizeUnit = "in" | "mm";
1821
@@ -23,13 +26,14 @@ interface FieldSpec {
2326 scale: number; // form value = SI value / scale
2427 step?: number;
2528 live?: boolean; // may be changed while a simulation runs
26- only?: Excitation; // shown only for this excitation
29+ show?: (p: PhysParams) => boolean; // shown only when this holds
2730 }
2831
2932 interface SelectSpec {
3033 key: keyof PhysParams;
3134 label: string;
3235 options: [string, string][];
36+ show?: (p: PhysParams) => boolean;
3337 }
3438
3539 interface Group {
@@ -37,9 +41,14 @@ interface Group {
3741 open: boolean;
3842 selects?: SelectSpec[]; // shown at the top of the group
3943 fields: FieldSpec[];
40- note?: string | "grid" | "damping";
44+ note?: string | "grid" | "damping" | "numerics" | "hits";
4145 }
4246
47+const push = (p: PhysParams) => p.excitation === "push";
48+const strike = (p: PhysParams) => p.excitation === "strike";
49+const series = (p: PhysParams) => strike(p) && p.hitSequence !== "single";
50+const graded = (p: PhysParams) => strike(p) && (p.hitSequence === "crescendo" || p.hitSequence === "decrescendo");
51+
4352 const groups: Group[] = [
4453 {
4554 title: "Sheet",
@@ -66,22 +75,29 @@ const groups: Group[] = [
6675 {
6776 title: "Excitation",
6877 open: true,
69- selects: [{ key: "excitation", label: "Type", options: [["push", "Push a corner and release"], ["strike", "Drumstick strike"]] }],
78+ selects: [
79+ { key: "excitation", label: "Type", options: [["strike", "Drumstick strike"], ["push", "Push a corner and release"]] },
80+ { key: "hitSequence", label: "Hits", options: [["single", "One hit"], ["steady", "Series, same speed"], ["crescendo", "Series, soft to loud"], ["decrescendo", "Series, loud to soft"]], show: strike },
81+ ],
7082 fields: [
71- { key: "pushAmount", label: "Push amount", unit: "mm", scale: 1e-3, step: 0.5, only: "push" },
72- { key: "fingerRadius", label: "Finger radius", unit: "mm", scale: 1e-3, step: 1, only: "push" },
73- { key: "rampTime", label: "Push ramp", unit: "s", scale: 1, step: 0.01, only: "push" },
74- { key: "holdMax", label: "Max settle time", unit: "s", scale: 1, step: 0.1, only: "push" },
75- { key: "releaseTime", label: "Release time", unit: "ms", scale: 1e-3, step: 1, only: "push" },
76- { key: "settleSigma0", label: "Settle damping σ₀", unit: "1/s", scale: 1, step: 10, only: "push" },
77- { key: "strikeFracX", label: "Strike x (fraction of width)", unit: "", scale: 1, step: 0.05, only: "strike" },
78- { key: "strikeFracY", label: "Strike y (fraction of height)", unit: "", scale: 1, step: 0.05, only: "strike" },
79- { key: "strikeVelocity", label: "Stick speed", unit: "m/s", scale: 1, step: 0.5, only: "strike" },
80- { key: "stickMass", label: "Stick mass", unit: "g", scale: 1e-3, step: 5, only: "strike" },
81- { key: "contactStiffness", label: "Contact stiffness K", unit: "N/mᵅ", scale: 1, step: 1e8, only: "strike" },
82- { key: "contactExponent", label: "Contact exponent α", unit: "", scale: 1, step: 0.1, only: "strike" },
83- { key: "contactRadius", label: "Tip radius", unit: "mm", scale: 1e-3, step: 1, only: "strike" },
83+ { key: "pushAmount", label: "Push amount", unit: "mm", scale: 1e-3, step: 0.5, show: push },
84+ { key: "fingerRadius", label: "Finger radius", unit: "mm", scale: 1e-3, step: 1, show: push },
85+ { key: "rampTime", label: "Push ramp", unit: "s", scale: 1, step: 0.01, show: push },
86+ { key: "holdMax", label: "Max settle time", unit: "s", scale: 1, step: 0.1, show: push },
87+ { key: "releaseTime", label: "Release time", unit: "ms", scale: 1e-3, step: 1, show: push },
88+ { key: "settleSigma0", label: "Settle damping σ₀", unit: "1/s", scale: 1, step: 10, show: push },
89+ { key: "strikeFracX", label: "Strike x (fraction of width)", unit: "", scale: 1, step: 0.05, show: strike },
90+ { key: "strikeFracY", label: "Strike y (fraction of height)", unit: "", scale: 1, step: 0.05, show: strike },
91+ { key: "strikeVelocity", label: "Stick speed", unit: "m/s", scale: 1, step: 0.5, show: strike },
92+ { key: "strikeVelocityMin", label: "Softest stick speed", unit: "m/s", scale: 1, step: 0.5, show: graded },
93+ { key: "hitCount", label: "Number of hits", unit: "", scale: 1, step: 1, show: series },
94+ { key: "hitInterval", label: "Time between hits", unit: "s", scale: 1, step: 0.1, show: series },
95+ { key: "stickMass", label: "Stick mass", unit: "g", scale: 1e-3, step: 5, show: strike },
96+ { key: "contactStiffness", label: "Contact stiffness K", unit: "N/mᵅ", scale: 1, step: 1e8, show: strike },
97+ { key: "contactExponent", label: "Contact exponent α", unit: "", scale: 1, step: 0.1, show: strike },
98+ { key: "contactRadius", label: "Tip radius", unit: "mm", scale: 1e-3, step: 1, show: strike },
8499 ],
100+ note: "hits",
85101 },
86102 {
87103 title: "Losses",
@@ -108,20 +124,22 @@ const groups: Group[] = [
108124 selects: [{ key: "rigidProjection", label: "Rigid-body projection while ringing", options: [["false", "Off"], ["true", "On (every 32 steps)"]] }],
109125 fields: [
110126 { key: "spacing", label: "Grid spacing", unit: "mm", scale: 1e-3, step: 0.5 },
127+ { key: "inPlaneMass", label: "In-plane mass factor", unit: "", scale: 1, step: 1 },
111128 { key: "fs", label: "Audio sample rate", unit: "Hz", scale: 1, step: 100 },
112129 { key: "fps", label: "Movie frame rate", unit: "fps", scale: 1, step: 10 },
113130 { key: "settleKE", label: "Settled below KE", unit: "J", scale: 1, step: 1e-9 },
114131 ],
115- note: "The excitation forces are self-equilibrated, so the projection only removes rounding drift of the six rigid modes. It is always applied while the sheet is held.",
132+ note: "numerics",
116133 },
117134 ];
118135
119136 const SIZE_SCALE: Record<SizeUnit, number> = { in: 0.0254, mm: 1e-3 };
120-const selectKeys: (keyof PhysParams)[] = ["shape", "bowShape", "excitation", "rigidProjection"];
137+const selectKeys: (keyof PhysParams)[] = ["shape", "bowShape", "excitation", "hitSequence", "rigidProjection"];
121138 const selectValues: Record<string, string[]> = {
122139 shape: ["rectangle", "ellipse"],
123140 bowShape: ["none", "width", "height", "dish"],
124141 excitation: ["push", "strike"],
142+ hitSequence: ["single", "steady", "crescendo", "decrescendo"],
125143 rigidProjection: ["false", "true"],
126144 };
127145
@@ -146,9 +164,12 @@ export class ParamForm {
146164 private inputs = new Map<keyof PhysParams, HTMLInputElement | HTMLSelectElement>();
147165 private sizeUnitSelect = document.createElement("select");
148166 private sizeUnitLabels: HTMLElement[] = [];
149- private conditional: { el: HTMLElement; only: Excitation }[] = [];
167+ private conditional: { el: HTMLElement; show: (p: PhysParams) => boolean }[] = [];
150168 private gridNote: HTMLElement | null = null;
169+ private hitsNote: HTMLElement | null = null;
151170 private dampingNote: HTMLElement | null = null;
171+ private numericsNote: HTMLElement | null = null;
172+ private omegaCache: { key: string; omega: number } | null = null;
152173 private liveKeys = new Set<keyof PhysParams>();
153174
154175 constructor(root: HTMLElement, params: PhysParams) {
@@ -159,10 +180,14 @@ export class ParamForm {
159180 const sum = document.createElement("summary");
160181 sum.textContent = g.title;
161182 det.appendChild(sum);
162- for (const s of g.selects ?? []) det.appendChild(this.selectRow(s));
183+ for (const s of g.selects ?? []) {
184+ const row = this.selectRow(s);
185+ if (s.show) this.conditional.push({ el: row, show: s.show });
186+ det.appendChild(row);
187+ }
163188 for (const f of g.fields) {
164189 const row = this.fieldRow(f);
165- if (f.only) this.conditional.push({ el: row, only: f.only });
190+ if (f.show) this.conditional.push({ el: row, show: f.show });
166191 if (f.live) this.liveKeys.add(f.key);
167192 det.appendChild(row);
168193 }
@@ -170,7 +195,9 @@ export class ParamForm {
170195 const note = document.createElement("p");
171196 note.className = "sub";
172197 if (g.note === "grid") this.gridNote = note;
198+ else if (g.note === "hits") this.hitsNote = note;
173199 else if (g.note === "damping") this.dampingNote = note;
200+ else if (g.note === "numerics") this.numericsNote = note;
174201 else note.textContent = g.note;
175202 det.appendChild(note);
176203 }
@@ -216,6 +243,7 @@ export class ParamForm {
216243 const v = parseFloat(inp.value);
217244 if (!isFinite(v)) return;
218245 (this.params as unknown as Record<string, number>)[f.key] = v * this.scale(f);
246+ this.updateVisibility();
219247 this.updateNotes();
220248 this.onChange(f.key, !!f.live);
221249 });
@@ -254,9 +282,9 @@ export class ParamForm {
254282 return f.unit === "size" ? SIZE_SCALE[this.sizeUnit] : f.scale;
255283 }
256284
257- /** Shows the fields of the selected excitation only. */
285+ /** Shows the fields that apply to the selected excitation and hit sequence only. */
258286 private updateVisibility(): void {
259- for (const c of this.conditional) c.el.hidden = c.only !== this.params.excitation;
287+ for (const c of this.conditional) c.el.hidden = !c.show(this.params);
260288 }
261289
262290 /** Rewrites every input from the current parameter values. */
@@ -306,8 +334,26 @@ export class ParamForm {
306334 return { grid, inside: countInside(grid) };
307335 }
308336
309- private updateNotes(): void {
337+ /** Rewrites the derived-quantity notes; call after a change made outside the form (the recorded duration). */
338+ updateNotes(): void {
310339 const p = this.params;
340+ if (this.hitsNote) {
341+ if (p.excitation !== "strike") {
342+ this.hitsNote.textContent = "";
343+ } else {
344+ const hits = hitSchedule(p);
345+ const last = hits[hits.length - 1];
346+ const decay = p.duration - last.t;
347+ const speeds = hits.length <= 8 ? hits.map((h) => h.v.toFixed(1)).join(", ") : `${hits[0].v.toFixed(1)} ... ${last.v.toFixed(1)}`;
348+ this.hitsNote.textContent =
349+ hits.length === 1
350+ ? `One hit at ${speeds} m/s; the recorded duration of ${p.duration} s starts at the impact.`
351+ : `${hits.length} hits at ${speeds} m/s, the last at ${last.t.toFixed(2)} s. ` +
352+ (decay > 0
353+ ? `The recorded duration of ${p.duration} s leaves ${decay.toFixed(2)} s of ring after the last hit.`
354+ : `The recorded duration of ${p.duration} s ends before the last hit; raise it to hear the whole series.`);
355+ }
356+ }
311357 if (this.gridNote) {
312358 const { grid, inside } = this.gridInfo();
313359 const m = makeMaterial(p);
@@ -320,6 +366,21 @@ export class ParamForm {
320366 const f = (hz: number) => t60(m, p, 2 * Math.PI * hz).toFixed(2);
321367 this.dampingNote.textContent = `T60: ${f(100)} s at 100 Hz, ${f(1000)} s at 1 kHz, ${f(5000)} s at 5 kHz.`;
322368 }
369+ if (this.numericsNote) {
370+ const { grid } = this.gridInfo();
371+ const m = makeMaterial(p);
372+ const s = Math.max(1, p.inPlaneMass);
373+ // The power iteration takes a moment, so its result is kept until the spacing, the material, or the factor changes.
374+ const key = [grid.hx, grid.hy, m.D, m.C, m.G, m.nu, m.rhoh, s].join("|");
375+ if (!this.omegaCache || this.omegaCache.key !== key) this.omegaCache = { key, omega: maxAngularFrequency(grid, m, s) };
376+ const stepsPerSample = Math.ceil(1 / (p.fs * ((CFL * 2) / this.omegaCache.omega)));
377+ const dt = 1 / (p.fs * stepsPerSample);
378+ const cp = Math.sqrt(m.C / m.rhoh / s);
379+ this.numericsNote.textContent =
380+ `Time step ${(dt * 1e6).toFixed(2)} µs (${stepsPerSample} per sample), set by the in-plane wave speed of ${cp.toFixed(0)} m/s. ` +
381+ `An in-plane mass factor above 1 slows the in-plane waves by its square root and lengthens the step by up to the same factor, at the cost of lowering the in-plane modes, which the bow couples to the sound. ` +
382+ `The excitation forces are self-equilibrated, so the projection while ringing only removes rounding drift of the six rigid modes; it is always applied while the sheet is held.`;
383+ }
323384 }
324385
325386 get excitation(): Excitation {
tests/gpu.test.tsmodified+11−1View file
@@ -34,6 +34,7 @@ interface Status {
3434 nx: number;
3535 ny: number;
3636 stick: { maxForce: number } | null;
37+ hits: { done: number; total: number } | null;
3738 }
3839
3940 const ARGS = ["--enable-unsafe-webgpu", "--enable-features=Vulkan", "--use-angle=vulkan", "--ignore-gpu-blocklist", "--no-sandbox"];
@@ -93,7 +94,7 @@ async function runToCompletion(query: string): Promise<Status> {
9394 }
9495
9596 test("a short push-and-release run completes with finite output", { timeout: 180_000 }, async () => {
96- const st = await runToCompletion("spacing=0.012&duration=0.2&rampTime=0.05&holdMax=0.1&rigidProjection=1");
97+ const st = await runToCompletion("spacing=0.012&duration=0.2&excitation=push&rampTime=0.05&holdMax=0.1&rigidProjection=1");
9798 console.log(` push: ${st.nx}x${st.ny}, dt ${(st.dt * 1e6).toFixed(2)} us, ${(st.stepsPerSec / 1e3).toFixed(0)} k steps/s, peak ${st.peakPressure.toFixed(3)} Pa, max |w| ${(st.wmax * 1e3).toFixed(2)} mm`);
9899 assert.equal(st.phase, "done");
99100 assert.ok(st.samples >= 0.2 * 48000 - 1, "recorded the requested duration");
@@ -109,3 +110,12 @@ test("a short drumstick strike on an elliptical sheet completes with finite outp
109110 assert.ok(st.peakPressure > 1e-3, "the sheet radiates");
110111 assert.ok(st.stick && st.stick.maxForce > 1, "the stick hit the sheet");
111112 });
113+
114+test("a short series of two strikes completes with finite output", { timeout: 180_000 }, async () => {
115+ const st = await runToCompletion("spacing=0.012&duration=0.2&hitSequence=crescendo&hitCount=2&hitInterval=0.1&strikeVelocityMin=1");
116+ console.log(` series: hits ${st.hits?.done} of ${st.hits?.total}, peak ${st.peakPressure.toFixed(3)} Pa, max contact force ${st.stick?.maxForce.toFixed(1)} N`);
117+ assert.equal(st.phase, "done");
118+ assert.deepEqual(st.hits, { done: 2, total: 2 });
119+ assert.ok(isFinite(st.wmax) && st.wmax < 0.05, "displacement stays bounded");
120+ assert.ok(st.stick && st.stick.maxForce > 1, "the stick hit the sheet");
121+});