5778036cymbalsim: scripted plate and cymbal designer with a modal linear model and a von Karman GPU modelJeremy Magland 1// Float64 reference implementation of the spatial discretization: the
2// forces of the geometrically nonlinear plate on the grid, the loss terms,
3// the discrete energy, and the largest frequency of the linearized system.
4// It mirrors the WGSL kernel of the full simulation and is used for the
5// time-step estimate, the tests, and the GPU-versus-CPU self-test; the
6// linear model assembles the same discretization as sparse matrices in
7// assemble.ts. It is not used in any time loop.
9import {
10 type Grid,
11 type Material,
12 type Plate,
13 type StandRates,
14 bendingStiffness,
15 cellThickness,
16 insideCell,
17 insideNode,
18 membraneStiffness,
19 nodeMass,
20 nodeWeight,
21 shearStiffness,
22 standRates,
23} from "../params.ts";
24import { type SigmaFit } from "../loss.ts";
26export interface Fields {
27 u: Float64Array;
28 v: Float64Array;
29 w: Float64Array;
30}
32export function fields(n: number): Fields {
33 return { u: new Float64Array(n), v: new Float64Array(n), w: new Float64Array(n) };
34}
36/** A plate of uniform thickness on the given grid (for the tests and the time-step estimate). */
37export function uniformPlate(g: Grid, mat: Material, h: number, w0?: Float64Array): Plate {
38 const hs = new Float64Array(g.n);
39 for (let j = 0; j < g.ny; j++) for (let i = 0; i < g.nx; i++) if (insideNode(g, i, j)) hs[j * g.nx + i] = h;
40 const freqs = new Float64Array([10, 40000]);
41 return { name: "uniform", Lx: (g.nx - 1) * g.hx, Ly: (g.ny - 1) * g.hy, grid: g, mat, h: hs, w0: w0 ?? new Float64Array(g.n), stand: null, loss: { freqs, t60: new Float64Array([1e9, 1e9]) } };
42}
44const G3 = 1 / Math.sqrt(3);
46/** Shape function derivatives of the bilinear cell at Gauss point gp (0..3), scaled by 2/h. */
47export function shapeDerivatives(gp: number, hx: number, hy: number, dNdx: Float64Array, dNdy: Float64Array): void {
48 const xi = gp & 1 ? G3 : -G3;
49 const eta = gp & 2 ? G3 : -G3;
50 dNdx[0] = (-(1 - eta) / 4) * (2 / hx);
51 dNdx[1] = ((1 - eta) / 4) * (2 / hx);
52 dNdx[2] = ((1 + eta) / 4) * (2 / hx);
53 dNdx[3] = (-(1 + eta) / 4) * (2 / hx);
54 dNdy[0] = (-(1 - xi) / 4) * (2 / hy);
55 dNdy[1] = (-(1 + xi) / 4) * (2 / hy);
56 dNdy[2] = ((1 + xi) / 4) * (2 / hy);
57 dNdy[3] = ((1 - xi) / 4) * (2 / hy);
58}
60/**
61 * Force of the pads at a node with transverse displacement w. Without
62 * lift-off the pads are a linear spring of rate kw. With lift-off each pad
63 * pushes only while compressed: the bottom pad by (preloadBottom - w), the
64 * top pad by (preloadTop + w); the constant part of the two, which balances
65 * the weight of the plate, is subtracted so that w = 0 is the rest state.
66 */
67export function padForce(sr: StandRates, k: number, w: number, liftoff: boolean): number {
68 if (!liftoff) return -sr.kw[k] * w;
69 const kb = sr.kwBottom[k];
70 const hasTop = sr.kw[k] > kb + 1e-300;
71 const bottom = kb * Math.max(sr.preloadBottom - w, 0);
72 const top = hasTop ? kb * Math.max(sr.preloadTop + w, 0) : 0;
73 const rest = kb * sr.preloadBottom - (hasTop ? kb * sr.preloadTop : 0);
74 return bottom - top - rest;
75}
77/**
78 * Undamped accelerations (out.u, out.v, out.w) of every node for the
79 * displacement fields f, measured from the stress-free shape pl.w0. Bending
80 * forces are the gradient of the discrete Kirchhoff energy with the local
81 * bending stiffness D(h) at the nodes and at the cells (natural free-edge
82 * conditions); membrane forces come from bilinear cells with 2x2 Gauss
83 * integration, the cell's mean thickness, and the Marguerre form of the von
84 * Karman strains; the pads of the stand add springs to ground over their
85 * annulus. Nodes outside the plate get zero.
86 */
87export function accelerations(pl: Plate, f: Fields, out: Fields, liftoff = false, sr: StandRates | null = standRates(pl)): void {
88 const g = pl.grid, m = pl.mat, w0 = pl.w0;
89 const { nx, ny, hx, hy } = g;
90 const { u, v, w } = f;
91 const hx2 = hx * hx, hy2 = hy * hy;
92 const nu = m.nu;
93 const id = (i: number, j: number) => j * nx + i;
94 const Dn = (i: number, j: number) => bendingStiffness(m, pl.h[id(i, j)]);
95 const Dc = (ci: number, cj: number) => bendingStiffness(m, cellThickness(pl, ci, cj));
97 const hasA = (i: number, j: number) => insideNode(g, i - 1, j) && insideNode(g, i, j) && insideNode(g, i + 1, j);
98 const hasB = (i: number, j: number) => insideNode(g, i, j - 1) && insideNode(g, i, j) && insideNode(g, i, j + 1);
99 const dxx = (i: number, j: number) => (w[id(i + 1, j)] - 2 * w[id(i, j)] + w[id(i - 1, j)]) / hx2;
100 const dyy = (i: number, j: number) => (w[id(i, j + 1)] - 2 * w[id(i, j)] + w[id(i, j - 1)]) / hy2;
101 const Ma = (i: number, j: number): number => {
102 if (!hasA(i, j)) return 0;
103 return nodeWeight(g, i, j) * Dn(i, j) * (dxx(i, j) + nu * (hasB(i, j) ? dyy(i, j) : 0));
104 };
105 const Mb = (i: number, j: number): number => {
106 if (!hasB(i, j)) return 0;
107 return nodeWeight(g, i, j) * Dn(i, j) * (dyy(i, j) + nu * (hasA(i, j) ? dxx(i, j) : 0));
108 };
109 const Cc = (ci: number, cj: number): number => {
110 if (!insideCell(g, ci, cj)) return 0;
111 return (Dc(ci, cj) * (w[id(ci + 1, cj + 1)] - w[id(ci + 1, cj)] - w[id(ci, cj + 1)] + w[id(ci, cj)])) / (hx * hy);
112 };
114 const q = (hx * hy) / 4;
115 const U = new Float64Array(4), V = new Float64Array(4), Wn = new Float64Array(4), W0 = new Float64Array(4);
116 const dNdx = new Float64Array(4), dNdy = new Float64Array(4);
118 const cellForce = (ci: number, cj: number, a: number, F: Float64Array): void => {
119 if (!insideCell(g, ci, cj)) return;
120 const hc = cellThickness(pl, ci, cj);
121 const C = membraneStiffness(m, hc), Gs = shearStiffness(m, hc);
122 const nn = [id(ci, cj), id(ci + 1, cj), id(ci + 1, cj + 1), id(ci, cj + 1)];
123 for (let b = 0; b < 4; b++) {
124 U[b] = u[nn[b]]; V[b] = v[nn[b]]; Wn[b] = w[nn[b]]; W0[b] = w0[nn[b]];
125 }
126 for (let gp = 0; gp < 4; gp++) {
127 shapeDerivatives(gp, hx, hy, dNdx, dNdy);
540490cForm the Marguerre strains from the deviation, removing a noise floor on curved designsJeremy Magland 128 let ux = 0, uy = 0, vx = 0, vy = 0, dwx = 0, dwy = 0, w0x = 0, w0y = 0;
5778036cymbalsim: scripted plate and cymbal designer with a modal linear model and a von Karman GPU modelJeremy Magland 129 for (let b = 0; b < 4; b++) {
130 ux += dNdx[b] * U[b]; uy += dNdy[b] * U[b];
131 vx += dNdx[b] * V[b]; vy += dNdy[b] * V[b];
540490cForm the Marguerre strains from the deviation, removing a noise floor on curved designsJeremy Magland 132 dwx += dNdx[b] * Wn[b]; dwy += dNdy[b] * Wn[b];
5778036cymbalsim: scripted plate and cymbal designer with a modal linear model and a von Karman GPU modelJeremy Magland 133 w0x += dNdx[b] * W0[b]; w0y += dNdy[b] * W0[b];
134 }
540490cForm the Marguerre strains from the deviation, removing a noise floor on curved designsJeremy Magland 135 const wx = w0x + dwx, wy = w0y + dwy;
136 // See step.wgsl: every term carries the deviation, so that the single
137 // precision of the GPU has nothing large to cancel.
138 const exx = ux + dwx * (w0x + 0.5 * dwx);
139 const eyy = vy + dwy * (w0y + 0.5 * dwy);
140 const gxy = uy + vx + w0x * dwy + w0y * dwx + dwx * dwy;
5778036cymbalsim: scripted plate and cymbal designer with a modal linear model and a von Karman GPU modelJeremy Magland 141 const Nxx = C * (exx + nu * eyy);
142 const Nyy = C * (eyy + nu * exx);
143 const Nxy = Gs * gxy;
144 F[0] -= q * (Nxx * dNdx[a] + Nxy * dNdy[a]);
145 F[1] -= q * (Nyy * dNdy[a] + Nxy * dNdx[a]);
146 F[2] -= q * ((Nxx * wx + Nxy * wy) * dNdx[a] + (Nxy * wx + Nyy * wy) * dNdy[a]);
147 }
148 };
150 const F = new Float64Array(3);
151 for (let j = 0; j < ny; j++) {
152 for (let i = 0; i < nx; i++) {
153 const k = id(i, j);
154 const mass = nodeMass(pl, i, j);
155 if (mass === 0) {
156 out.u[k] = 0; out.v[k] = 0; out.w[k] = 0;
157 continue;
158 }
159 const bend =
160 (Ma(i - 1, j) - 2 * Ma(i, j) + Ma(i + 1, j)) / hx2 +
161 (Mb(i, j - 1) - 2 * Mb(i, j) + Mb(i, j + 1)) / hy2 +
162 (2 * (1 - nu) * (Cc(i, j) - Cc(i - 1, j) - Cc(i, j - 1) + Cc(i - 1, j - 1))) / (hx * hy);
163 F[0] = 0; F[1] = 0; F[2] = -hx * hy * bend;
164 if (sr) {
165 F[0] -= sr.ku[k] * u[k]; F[1] -= sr.ku[k] * v[k];
166 F[2] += padForce(sr, k, w[k], liftoff);
167 }
168 cellForce(i, j, 0, F);
169 cellForce(i - 1, j, 1, F);
170 cellForce(i - 1, j - 1, 2, F);
171 cellForce(i, j - 1, 3, F);
172 out.u[k] = F[0] / mass;
173 out.v[k] = F[1] / mass;
174 out.w[k] = F[2] / mass;
175 }
176 }
177}
179/** The loss terms of the full simulation: the fitted sigma model on the transverse motion, an in-plane rate, and the pads' dashpots. */
180export interface FullLoss {
181 sigma0: number; // [1/s], applied as 2 sigma0 on all components
182 sigma1: number; // [m^2/s], applied as 2 sigma1 times the Laplacian of the transverse velocity
183 sigmaIn: number; // [1/s], extra rate on the in-plane motion
184 /** Dashpot rate of the pads per unit of pad stiffness [s]: c_k = padDamping * kw_k. */
185 padDamping: number;
186}
188/** The loss terms of the full simulation from the fitted sigma model and the stand. */
189export function fullLoss(fit: SigmaFit, sr: StandRates | null, inPlaneRate: number): FullLoss {
190 const padDamping = sr && sr.bounceHz > 0 ? sr.eta / (2 * Math.PI * sr.bounceHz) : 0;
191 return { sigma0: fit.sigma0, sigma1: fit.sigma1, sigmaIn: inPlaneRate, padDamping };
192}
194/**
195 * Adds the loss terms to the accelerations `out` for the velocities `vel`:
196 * -2 sigma0 v on every component, +2 sigma1 (graph Laplacian of w_t) on the
197 * transverse component, -2 sigmaIn v on the in-plane components, and the
198 * dashpots of the pads over their annulus.
199 */
200export function addDamping(pl: Plate, loss: FullLoss, vel: Fields, out: Fields, sr: StandRates | null = standRates(pl)): void {
201 const g = pl.grid;
202 const { sigma0, sigma1, sigmaIn } = loss;
203 const { nx, ny, hx, hy } = g;
204 const id = (i: number, j: number) => j * nx + i;
205 for (let j = 0; j < ny; j++) {
206 for (let i = 0; i < nx; i++) {
207 const k = id(i, j);
208 const mass = nodeMass(pl, i, j);
209 if (mass === 0) continue;
210 const c = vel.w[k];
211 let lap = 0;
212 if (insideNode(g, i - 1, j)) lap += (vel.w[id(i - 1, j)] - c) / (hx * hx);
213 if (insideNode(g, i + 1, j)) lap += (vel.w[id(i + 1, j)] - c) / (hx * hx);
214 if (insideNode(g, i, j - 1)) lap += (vel.w[id(i, j - 1)] - c) / (hy * hy);
215 if (insideNode(g, i, j + 1)) lap += (vel.w[id(i, j + 1)] - c) / (hy * hy);
216 const cw = sr ? (loss.padDamping * sr.kw[k]) / mass : 0;
217 const cu = sr ? (loss.padDamping * sr.ku[k]) / mass : 0;
218 out.u[k] -= (2 * (sigma0 + sigmaIn) + cu) * vel.u[k];
219 out.v[k] -= (2 * (sigma0 + sigmaIn) + cu) * vel.v[k];
220 out.w[k] += -(2 * sigma0 + cw) * c + 2 * sigma1 * lap;
221 }
222 }
223}
225/** Discrete potential energy (bending + membrane + pads, linear pads) whose negative gradient is the force in `accelerations`. */
226export function potentialEnergy(pl: Plate, f: Fields, sr: StandRates | null = standRates(pl)): number {
227 const g = pl.grid, m = pl.mat, w0 = pl.w0;
228 const { nx, ny, hx, hy } = g;
229 const { u, v, w } = f;
230 const nu = m.nu;
231 const id = (i: number, j: number) => j * nx + i;
232 let Ub = 0;
233 for (let j = 0; j < ny; j++) {
234 for (let i = 0; i < nx; i++) {
235 let a = 0, b = 0;
236 if (insideNode(g, i - 1, j) && insideNode(g, i, j) && insideNode(g, i + 1, j)) a = (w[id(i + 1, j)] - 2 * w[id(i, j)] + w[id(i - 1, j)]) / (hx * hx);
237 if (insideNode(g, i, j - 1) && insideNode(g, i, j) && insideNode(g, i, j + 1)) b = (w[id(i, j + 1)] - 2 * w[id(i, j)] + w[id(i, j - 1)]) / (hy * hy);
238 Ub += nodeWeight(g, i, j) * bendingStiffness(m, pl.h[id(i, j)]) * (a * a + b * b + 2 * nu * a * b);
239 }
240 }
241 for (let cj = 0; cj < ny - 1; cj++) {
242 for (let ci = 0; ci < nx - 1; ci++) {
243 if (!insideCell(g, ci, cj)) continue;
244 const c = (w[id(ci + 1, cj + 1)] - w[id(ci + 1, cj)] - w[id(ci, cj + 1)] + w[id(ci, cj)]) / (hx * hy);
245 Ub += 2 * (1 - nu) * bendingStiffness(m, cellThickness(pl, ci, cj)) * c * c;
246 }
247 }
248 Ub *= 0.5 * hx * hy;
250 let Um = 0;
251 const q = (hx * hy) / 4;
252 const dNdx = new Float64Array(4), dNdy = new Float64Array(4);
253 for (let cj = 0; cj < ny - 1; cj++) {
254 for (let ci = 0; ci < nx - 1; ci++) {
255 if (!insideCell(g, ci, cj)) continue;
256 const hc = cellThickness(pl, ci, cj);
257 const C = membraneStiffness(m, hc), Gs = shearStiffness(m, hc);
258 const nn = [id(ci, cj), id(ci + 1, cj), id(ci + 1, cj + 1), id(ci, cj + 1)];
259 for (let gp = 0; gp < 4; gp++) {
260 shapeDerivatives(gp, hx, hy, dNdx, dNdy);
540490cForm the Marguerre strains from the deviation, removing a noise floor on curved designsJeremy Magland 261 let ux = 0, uy = 0, vx = 0, vy = 0, dwx = 0, dwy = 0, w0x = 0, w0y = 0;
5778036cymbalsim: scripted plate and cymbal designer with a modal linear model and a von Karman GPU modelJeremy Magland 262 for (let b = 0; b < 4; b++) {
263 const z0 = w0[nn[b]];
264 ux += dNdx[b] * u[nn[b]]; uy += dNdy[b] * u[nn[b]];
265 vx += dNdx[b] * v[nn[b]]; vy += dNdy[b] * v[nn[b]];
540490cForm the Marguerre strains from the deviation, removing a noise floor on curved designsJeremy Magland 266 dwx += dNdx[b] * w[nn[b]]; dwy += dNdy[b] * w[nn[b]];
5778036cymbalsim: scripted plate and cymbal designer with a modal linear model and a von Karman GPU modelJeremy Magland 267 w0x += dNdx[b] * z0; w0y += dNdy[b] * z0;
268 }
540490cForm the Marguerre strains from the deviation, removing a noise floor on curved designsJeremy Magland 269 const exx = ux + dwx * (w0x + 0.5 * dwx);
270 const eyy = vy + dwy * (w0y + 0.5 * dwy);
271 const gxy = uy + vx + w0x * dwy + w0y * dwx + dwx * dwy;
5778036cymbalsim: scripted plate and cymbal designer with a modal linear model and a von Karman GPU modelJeremy Magland 272 Um += q * 0.5 * (C * (exx * exx + eyy * eyy + 2 * nu * exx * eyy) + Gs * gxy * gxy);
273 }
274 }
275 }
277 let Us = 0;
278 if (sr) for (const k of sr.nodes) Us += 0.5 * (sr.ku[k] * (u[k] * u[k] + v[k] * v[k]) + sr.kw[k] * w[k] * w[k]);
279 return Ub + Um + Us;
280}
282/**
283 * Largest angular frequency of the linearized discrete system, by power
284 * iteration on a small full grid with the same spacing and a uniform
285 * thickness h (the largest of the plate, which bounds the bending frequency;
286 * the in-plane wave speed does not depend on the thickness). The extreme
287 * mode is a grid-scale mode, so neither the plate size nor its outline
288 * matters.
289 */
290export function maxAngularFrequency(g: Grid, m: Material, h: number, nSmall = 12, iters = 300): number {
291 const gs: Grid = { nx: nSmall, ny: nSmall, hx: g.hx, hy: g.hy, n: nSmall * nSmall };
292 const pl = uniformPlate(gs, m, h);
293 const n = gs.n;
294 const x = fields(n), y = fields(n);
295 let seed = 12345;
296 const rnd = () => {
297 seed = (seed * 1664525 + 1013904223) >>> 0;
298 return seed / 4294967296 - 0.5;
299 };
300 for (let k = 0; k < n; k++) {
301 x.u[k] = rnd(); x.v[k] = rnd(); x.w[k] = rnd();
302 }
303 const eps = 1e-9;
304 let lambda = 0;
305 for (let it = 0; it < iters; it++) {
306 let norm = 0;
307 for (let k = 0; k < n; k++) norm += x.u[k] ** 2 + x.v[k] ** 2 + x.w[k] ** 2;
308 norm = Math.sqrt(norm);
309 for (let k = 0; k < n; k++) {
310 x.u[k] *= eps / norm; x.v[k] *= eps / norm; x.w[k] *= eps / norm;
311 }
312 accelerations(pl, x, y, false, null);
313 let dot = 0;
314 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];
315 lambda = -dot / (eps * eps);
316 for (let k = 0; k < n; k++) {
317 x.u[k] = -y.u[k]; x.v[k] = -y.v[k]; x.w[k] = -y.w[k];
318 }
319 }
320 return Math.sqrt(Math.max(lambda, 0));
321}