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);
128 let ux = 0, uy = 0, vx = 0, vy = 0, wx = 0, wy = 0, w0x = 0, w0y = 0;
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];
132 wx += dNdx[b] * (Wn[b] + W0[b]); wy += dNdy[b] * (Wn[b] + W0[b]);
133 w0x += dNdx[b] * W0[b]; w0y += dNdy[b] * W0[b];
134 }
135 const exx = ux + 0.5 * (wx * wx - w0x * w0x);
136 const eyy = vy + 0.5 * (wy * wy - w0y * w0y);
137 const gxy = uy + vx + wx * wy - w0x * w0y;
138 const Nxx = C * (exx + nu * eyy);
139 const Nyy = C * (eyy + nu * exx);
140 const Nxy = Gs * gxy;
141 F[0] -= q * (Nxx * dNdx[a] + Nxy * dNdy[a]);
142 F[1] -= q * (Nyy * dNdy[a] + Nxy * dNdx[a]);
143 F[2] -= q * ((Nxx * wx + Nxy * wy) * dNdx[a] + (Nxy * wx + Nyy * wy) * dNdy[a]);
144 }
145 };
147 const F = new Float64Array(3);
148 for (let j = 0; j < ny; j++) {
149 for (let i = 0; i < nx; i++) {
150 const k = id(i, j);
151 const mass = nodeMass(pl, i, j);
152 if (mass === 0) {
153 out.u[k] = 0; out.v[k] = 0; out.w[k] = 0;
154 continue;
155 }
156 const bend =
157 (Ma(i - 1, j) - 2 * Ma(i, j) + Ma(i + 1, j)) / hx2 +
158 (Mb(i, j - 1) - 2 * Mb(i, j) + Mb(i, j + 1)) / hy2 +
159 (2 * (1 - nu) * (Cc(i, j) - Cc(i - 1, j) - Cc(i, j - 1) + Cc(i - 1, j - 1))) / (hx * hy);
160 F[0] = 0; F[1] = 0; F[2] = -hx * hy * bend;
161 if (sr) {
162 F[0] -= sr.ku[k] * u[k]; F[1] -= sr.ku[k] * v[k];
163 F[2] += padForce(sr, k, w[k], liftoff);
164 }
165 cellForce(i, j, 0, F);
166 cellForce(i - 1, j, 1, F);
167 cellForce(i - 1, j - 1, 2, F);
168 cellForce(i, j - 1, 3, F);
169 out.u[k] = F[0] / mass;
170 out.v[k] = F[1] / mass;
171 out.w[k] = F[2] / mass;
172 }
173 }
174}
176/** The loss terms of the full simulation: the fitted sigma model on the transverse motion, an in-plane rate, and the pads' dashpots. */
177export interface FullLoss {
178 sigma0: number; // [1/s], applied as 2 sigma0 on all components
179 sigma1: number; // [m^2/s], applied as 2 sigma1 times the Laplacian of the transverse velocity
180 sigmaIn: number; // [1/s], extra rate on the in-plane motion
181 /** Dashpot rate of the pads per unit of pad stiffness [s]: c_k = padDamping * kw_k. */
182 padDamping: number;
183}
185/** The loss terms of the full simulation from the fitted sigma model and the stand. */
186export function fullLoss(fit: SigmaFit, sr: StandRates | null, inPlaneRate: number): FullLoss {
187 const padDamping = sr && sr.bounceHz > 0 ? sr.eta / (2 * Math.PI * sr.bounceHz) : 0;
188 return { sigma0: fit.sigma0, sigma1: fit.sigma1, sigmaIn: inPlaneRate, padDamping };
189}
191/**
192 * Adds the loss terms to the accelerations `out` for the velocities `vel`:
193 * -2 sigma0 v on every component, +2 sigma1 (graph Laplacian of w_t) on the
194 * transverse component, -2 sigmaIn v on the in-plane components, and the
195 * dashpots of the pads over their annulus.
196 */
197export function addDamping(pl: Plate, loss: FullLoss, vel: Fields, out: Fields, sr: StandRates | null = standRates(pl)): void {
198 const g = pl.grid;
199 const { sigma0, sigma1, sigmaIn } = loss;
200 const { nx, ny, hx, hy } = g;
201 const id = (i: number, j: number) => j * nx + i;
202 for (let j = 0; j < ny; j++) {
203 for (let i = 0; i < nx; i++) {
204 const k = id(i, j);
205 const mass = nodeMass(pl, i, j);
206 if (mass === 0) continue;
207 const c = vel.w[k];
208 let lap = 0;
209 if (insideNode(g, i - 1, j)) lap += (vel.w[id(i - 1, j)] - c) / (hx * hx);
210 if (insideNode(g, i + 1, j)) lap += (vel.w[id(i + 1, j)] - c) / (hx * hx);
211 if (insideNode(g, i, j - 1)) lap += (vel.w[id(i, j - 1)] - c) / (hy * hy);
212 if (insideNode(g, i, j + 1)) lap += (vel.w[id(i, j + 1)] - c) / (hy * hy);
213 const cw = sr ? (loss.padDamping * sr.kw[k]) / mass : 0;
214 const cu = sr ? (loss.padDamping * sr.ku[k]) / mass : 0;
215 out.u[k] -= (2 * (sigma0 + sigmaIn) + cu) * vel.u[k];
216 out.v[k] -= (2 * (sigma0 + sigmaIn) + cu) * vel.v[k];
217 out.w[k] += -(2 * sigma0 + cw) * c + 2 * sigma1 * lap;
218 }
219 }
220}
222/** Discrete potential energy (bending + membrane + pads, linear pads) whose negative gradient is the force in `accelerations`. */
223export function potentialEnergy(pl: Plate, f: Fields, sr: StandRates | null = standRates(pl)): number {
224 const g = pl.grid, m = pl.mat, w0 = pl.w0;
225 const { nx, ny, hx, hy } = g;
226 const { u, v, w } = f;
227 const nu = m.nu;
228 const id = (i: number, j: number) => j * nx + i;
229 let Ub = 0;
230 for (let j = 0; j < ny; j++) {
231 for (let i = 0; i < nx; i++) {
232 let a = 0, b = 0;
233 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);
234 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);
235 Ub += nodeWeight(g, i, j) * bendingStiffness(m, pl.h[id(i, j)]) * (a * a + b * b + 2 * nu * a * b);
236 }
237 }
238 for (let cj = 0; cj < ny - 1; cj++) {
239 for (let ci = 0; ci < nx - 1; ci++) {
240 if (!insideCell(g, ci, cj)) continue;
241 const c = (w[id(ci + 1, cj + 1)] - w[id(ci + 1, cj)] - w[id(ci, cj + 1)] + w[id(ci, cj)]) / (hx * hy);
242 Ub += 2 * (1 - nu) * bendingStiffness(m, cellThickness(pl, ci, cj)) * c * c;
243 }
244 }
245 Ub *= 0.5 * hx * hy;
247 let Um = 0;
248 const q = (hx * hy) / 4;
249 const dNdx = new Float64Array(4), dNdy = new Float64Array(4);
250 for (let cj = 0; cj < ny - 1; cj++) {
251 for (let ci = 0; ci < nx - 1; ci++) {
252 if (!insideCell(g, ci, cj)) continue;
253 const hc = cellThickness(pl, ci, cj);
254 const C = membraneStiffness(m, hc), Gs = shearStiffness(m, hc);
255 const nn = [id(ci, cj), id(ci + 1, cj), id(ci + 1, cj + 1), id(ci, cj + 1)];
256 for (let gp = 0; gp < 4; gp++) {
257 shapeDerivatives(gp, hx, hy, dNdx, dNdy);
258 let ux = 0, uy = 0, vx = 0, vy = 0, wx = 0, wy = 0, w0x = 0, w0y = 0;
259 for (let b = 0; b < 4; b++) {
260 const z0 = w0[nn[b]];
261 ux += dNdx[b] * u[nn[b]]; uy += dNdy[b] * u[nn[b]];
262 vx += dNdx[b] * v[nn[b]]; vy += dNdy[b] * v[nn[b]];
263 wx += dNdx[b] * (w[nn[b]] + z0); wy += dNdy[b] * (w[nn[b]] + z0);
264 w0x += dNdx[b] * z0; w0y += dNdy[b] * z0;
265 }
266 const exx = ux + 0.5 * (wx * wx - w0x * w0x);
267 const eyy = vy + 0.5 * (wy * wy - w0y * w0y);
268 const gxy = uy + vx + wx * wy - w0x * w0y;
269 Um += q * 0.5 * (C * (exx * exx + eyy * eyy + 2 * nu * exx * eyy) + Gs * gxy * gxy);
270 }
271 }
272 }
274 let Us = 0;
275 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]);
276 return Ub + Um + Us;
277}
279/**
280 * Largest angular frequency of the linearized discrete system, by power
281 * iteration on a small full grid with the same spacing and a uniform
282 * thickness h (the largest of the plate, which bounds the bending frequency;
283 * the in-plane wave speed does not depend on the thickness). The extreme
284 * mode is a grid-scale mode, so neither the plate size nor its outline
285 * matters.
286 */
287export function maxAngularFrequency(g: Grid, m: Material, h: number, nSmall = 12, iters = 300): number {
288 const gs: Grid = { nx: nSmall, ny: nSmall, hx: g.hx, hy: g.hy, n: nSmall * nSmall };
289 const pl = uniformPlate(gs, m, h);
290 const n = gs.n;
291 const x = fields(n), y = fields(n);
292 let seed = 12345;
293 const rnd = () => {
294 seed = (seed * 1664525 + 1013904223) >>> 0;
295 return seed / 4294967296 - 0.5;
296 };
297 for (let k = 0; k < n; k++) {
298 x.u[k] = rnd(); x.v[k] = rnd(); x.w[k] = rnd();
299 }
300 const eps = 1e-9;
301 let lambda = 0;
302 for (let it = 0; it < iters; it++) {
303 let norm = 0;
304 for (let k = 0; k < n; k++) norm += x.u[k] ** 2 + x.v[k] ** 2 + x.w[k] ** 2;
305 norm = Math.sqrt(norm);
306 for (let k = 0; k < n; k++) {
307 x.u[k] *= eps / norm; x.v[k] *= eps / norm; x.w[k] *= eps / norm;
308 }
309 accelerations(pl, x, y, false, null);
310 let dot = 0;
311 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];
312 lambda = -dot / (eps * eps);
313 for (let k = 0; k < n; k++) {
314 x.u[k] = -y.u[k]; x.v[k] = -y.v[k]; x.w[k] = -y.w[k];
315 }
316 }
317 return Math.sqrt(Math.max(lambda, 0));
318}