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,
2a1a6daDerive the losses from physical constants instead of a prescribed decay curveJeremy Magland 13 ROOM_TEMPERATURE,
5778036cymbalsim: scripted plate and cymbal designer with a modal linear model and a von Karman GPU modelJeremy Magland 14 bendingStiffness,
15 cellThickness,
16 insideCell,
17 insideNode,
18 membraneStiffness,
19 nodeMass,
20 nodeWeight,
21 shearStiffness,
22} from "../params.ts";
23import { type SigmaFit } from "../loss.ts";
25export interface Fields {
26 u: Float64Array;
27 v: Float64Array;
28 w: Float64Array;
29}
31export function fields(n: number): Fields {
32 return { u: new Float64Array(n), v: new Float64Array(n), w: new Float64Array(n) };
33}
35/** A plate of uniform thickness on the given grid (for the tests and the time-step estimate). */
36export function uniformPlate(g: Grid, mat: Material, h: number, w0?: Float64Array): Plate {
37 const hs = new Float64Array(g.n);
38 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;
2481a63Remove the standJeremy Magland 39 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), temperature: ROOM_TEMPERATURE };
5778036cymbalsim: scripted plate and cymbal designer with a modal linear model and a von Karman GPU modelJeremy Magland 40}
42const G3 = 1 / Math.sqrt(3);
44/** Shape function derivatives of the bilinear cell at Gauss point gp (0..3), scaled by 2/h. */
45export function shapeDerivatives(gp: number, hx: number, hy: number, dNdx: Float64Array, dNdy: Float64Array): void {
46 const xi = gp & 1 ? G3 : -G3;
47 const eta = gp & 2 ? G3 : -G3;
48 dNdx[0] = (-(1 - eta) / 4) * (2 / hx);
49 dNdx[1] = ((1 - eta) / 4) * (2 / hx);
50 dNdx[2] = ((1 + eta) / 4) * (2 / hx);
51 dNdx[3] = (-(1 + eta) / 4) * (2 / hx);
52 dNdy[0] = (-(1 - xi) / 4) * (2 / hy);
53 dNdy[1] = (-(1 + xi) / 4) * (2 / hy);
54 dNdy[2] = ((1 + xi) / 4) * (2 / hy);
55 dNdy[3] = ((1 - xi) / 4) * (2 / hy);
56}
58/**
59 * Undamped accelerations (out.u, out.v, out.w) of every node for the
60 * displacement fields f, measured from the stress-free shape pl.w0. Bending
61 * forces are the gradient of the discrete Kirchhoff energy with the local
62 * bending stiffness D(h) at the nodes and at the cells (natural free-edge
63 * conditions); membrane forces come from bilinear cells with 2x2 Gauss
64 * integration, the cell's mean thickness, and the Marguerre form of the von
5778036cymbalsim: scripted plate and cymbal designer with a modal linear model and a von Karman GPU modelJeremy Magland 66 */
2481a63Remove the standJeremy Magland 67export function accelerations(pl: Plate, f: Fields, out: Fields): void {
5778036cymbalsim: scripted plate and cymbal designer with a modal linear model and a von Karman GPU modelJeremy Magland 68 const g = pl.grid, m = pl.mat, w0 = pl.w0;
69 const { nx, ny, hx, hy } = g;
70 const { u, v, w } = f;
71 const hx2 = hx * hx, hy2 = hy * hy;
72 const nu = m.nu;
73 const id = (i: number, j: number) => j * nx + i;
74 const Dn = (i: number, j: number) => bendingStiffness(m, pl.h[id(i, j)]);
75 const Dc = (ci: number, cj: number) => bendingStiffness(m, cellThickness(pl, ci, cj));
77 const hasA = (i: number, j: number) => insideNode(g, i - 1, j) && insideNode(g, i, j) && insideNode(g, i + 1, j);
78 const hasB = (i: number, j: number) => insideNode(g, i, j - 1) && insideNode(g, i, j) && insideNode(g, i, j + 1);
79 const dxx = (i: number, j: number) => (w[id(i + 1, j)] - 2 * w[id(i, j)] + w[id(i - 1, j)]) / hx2;
80 const dyy = (i: number, j: number) => (w[id(i, j + 1)] - 2 * w[id(i, j)] + w[id(i, j - 1)]) / hy2;
81 const Ma = (i: number, j: number): number => {
82 if (!hasA(i, j)) return 0;
83 return nodeWeight(g, i, j) * Dn(i, j) * (dxx(i, j) + nu * (hasB(i, j) ? dyy(i, j) : 0));
84 };
85 const Mb = (i: number, j: number): number => {
86 if (!hasB(i, j)) return 0;
87 return nodeWeight(g, i, j) * Dn(i, j) * (dyy(i, j) + nu * (hasA(i, j) ? dxx(i, j) : 0));
88 };
89 const Cc = (ci: number, cj: number): number => {
90 if (!insideCell(g, ci, cj)) return 0;
91 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);
92 };
94 const q = (hx * hy) / 4;
95 const U = new Float64Array(4), V = new Float64Array(4), Wn = new Float64Array(4), W0 = new Float64Array(4);
96 const dNdx = new Float64Array(4), dNdy = new Float64Array(4);
98 const cellForce = (ci: number, cj: number, a: number, F: Float64Array): void => {
99 if (!insideCell(g, ci, cj)) return;
100 const hc = cellThickness(pl, ci, cj);
101 const C = membraneStiffness(m, hc), Gs = shearStiffness(m, hc);
102 const nn = [id(ci, cj), id(ci + 1, cj), id(ci + 1, cj + 1), id(ci, cj + 1)];
103 for (let b = 0; b < 4; b++) {
104 U[b] = u[nn[b]]; V[b] = v[nn[b]]; Wn[b] = w[nn[b]]; W0[b] = w0[nn[b]];
105 }
106 for (let gp = 0; gp < 4; gp++) {
107 shapeDerivatives(gp, hx, hy, dNdx, dNdy);
540490cForm the Marguerre strains from the deviation, removing a noise floor on curved designsJeremy Magland 108 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 109 for (let b = 0; b < 4; b++) {
110 ux += dNdx[b] * U[b]; uy += dNdy[b] * U[b];
111 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 112 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 113 w0x += dNdx[b] * W0[b]; w0y += dNdy[b] * W0[b];
114 }
540490cForm the Marguerre strains from the deviation, removing a noise floor on curved designsJeremy Magland 115 const wx = w0x + dwx, wy = w0y + dwy;
116 // See step.wgsl: every term carries the deviation, so that the single
117 // precision of the GPU has nothing large to cancel.
118 const exx = ux + dwx * (w0x + 0.5 * dwx);
119 const eyy = vy + dwy * (w0y + 0.5 * dwy);
120 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 121 const Nxx = C * (exx + nu * eyy);
122 const Nyy = C * (eyy + nu * exx);
123 const Nxy = Gs * gxy;
124 F[0] -= q * (Nxx * dNdx[a] + Nxy * dNdy[a]);
125 F[1] -= q * (Nyy * dNdy[a] + Nxy * dNdx[a]);
126 F[2] -= q * ((Nxx * wx + Nxy * wy) * dNdx[a] + (Nxy * wx + Nyy * wy) * dNdy[a]);
127 }
128 };
130 const F = new Float64Array(3);
131 for (let j = 0; j < ny; j++) {
132 for (let i = 0; i < nx; i++) {
133 const k = id(i, j);
134 const mass = nodeMass(pl, i, j);
135 if (mass === 0) {
136 out.u[k] = 0; out.v[k] = 0; out.w[k] = 0;
137 continue;
138 }
139 const bend =
140 (Ma(i - 1, j) - 2 * Ma(i, j) + Ma(i + 1, j)) / hx2 +
141 (Mb(i, j - 1) - 2 * Mb(i, j) + Mb(i, j + 1)) / hy2 +
142 (2 * (1 - nu) * (Cc(i, j) - Cc(i - 1, j) - Cc(i, j - 1) + Cc(i - 1, j - 1))) / (hx * hy);
143 F[0] = 0; F[1] = 0; F[2] = -hx * hy * bend;
144 cellForce(i, j, 0, F);
145 cellForce(i - 1, j, 1, F);
146 cellForce(i - 1, j - 1, 2, F);
147 cellForce(i, j - 1, 3, F);
148 out.u[k] = F[0] / mass;
149 out.v[k] = F[1] / mass;
150 out.w[k] = F[2] / mass;
151 }
152 }
153}
2481a63Remove the standJeremy Magland 155/** The loss terms of the full simulation: the fitted sigma model on the transverse motion and an in-plane rate. */
5778036cymbalsim: scripted plate and cymbal designer with a modal linear model and a von Karman GPU modelJeremy Magland 156export interface FullLoss {
157 sigma0: number; // [1/s], applied as 2 sigma0 on all components
158 sigma1: number; // [m^2/s], applied as 2 sigma1 times the Laplacian of the transverse velocity
159 sigmaIn: number; // [1/s], extra rate on the in-plane motion
160}
18b1b86Model the felt as a Maxwell chain so that the full model has the loss factor of the linear oneJeremy Magland 162/** The loss terms of the full simulation from the fitted sigma model. */
163export function fullLoss(fit: SigmaFit, inPlaneRate: number): FullLoss {
164 return { sigma0: fit.sigma0, sigma1: fit.sigma1, sigmaIn: inPlaneRate };
5778036cymbalsim: scripted plate and cymbal designer with a modal linear model and a von Karman GPU modelJeremy Magland 165}
167/**
168 * Adds the loss terms to the accelerations `out` for the velocities `vel`:
169 * -2 sigma0 v on every component, +2 sigma1 (graph Laplacian of w_t) on the
18b1b86Model the felt as a Maxwell chain so that the full model has the loss factor of the linear oneJeremy Magland 170 * transverse component, and -2 sigmaIn v on the in-plane components.
5778036cymbalsim: scripted plate and cymbal designer with a modal linear model and a von Karman GPU modelJeremy Magland 171 */
18b1b86Model the felt as a Maxwell chain so that the full model has the loss factor of the linear oneJeremy Magland 172export function addDamping(pl: Plate, loss: FullLoss, vel: Fields, out: Fields): void {
5778036cymbalsim: scripted plate and cymbal designer with a modal linear model and a von Karman GPU modelJeremy Magland 173 const g = pl.grid;
174 const { sigma0, sigma1, sigmaIn } = loss;
175 const { nx, ny, hx, hy } = g;
176 const id = (i: number, j: number) => j * nx + i;
177 for (let j = 0; j < ny; j++) {
178 for (let i = 0; i < nx; i++) {
179 const k = id(i, j);
180 const mass = nodeMass(pl, i, j);
181 if (mass === 0) continue;
182 const c = vel.w[k];
183 let lap = 0;
184 if (insideNode(g, i - 1, j)) lap += (vel.w[id(i - 1, j)] - c) / (hx * hx);
185 if (insideNode(g, i + 1, j)) lap += (vel.w[id(i + 1, j)] - c) / (hx * hx);
186 if (insideNode(g, i, j - 1)) lap += (vel.w[id(i, j - 1)] - c) / (hy * hy);
187 if (insideNode(g, i, j + 1)) lap += (vel.w[id(i, j + 1)] - c) / (hy * hy);
18b1b86Model the felt as a Maxwell chain so that the full model has the loss factor of the linear oneJeremy Magland 188 out.u[k] -= 2 * (sigma0 + sigmaIn) * vel.u[k];
189 out.v[k] -= 2 * (sigma0 + sigmaIn) * vel.v[k];
190 out.w[k] += -2 * sigma0 * c + 2 * sigma1 * lap;
5778036cymbalsim: scripted plate and cymbal designer with a modal linear model and a von Karman GPU modelJeremy Magland 191 }
192 }
193}
2481a63Remove the standJeremy Magland 195/** Discrete potential energy (bending + membrane) whose negative gradient is the force in `accelerations`. */
196export function potentialEnergy(pl: Plate, f: Fields): number {
5778036cymbalsim: scripted plate and cymbal designer with a modal linear model and a von Karman GPU modelJeremy Magland 197 const g = pl.grid, m = pl.mat, w0 = pl.w0;
198 const { nx, ny, hx, hy } = g;
199 const { u, v, w } = f;
200 const nu = m.nu;
201 const id = (i: number, j: number) => j * nx + i;
202 let Ub = 0;
203 for (let j = 0; j < ny; j++) {
204 for (let i = 0; i < nx; i++) {
205 let a = 0, b = 0;
206 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);
207 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);
208 Ub += nodeWeight(g, i, j) * bendingStiffness(m, pl.h[id(i, j)]) * (a * a + b * b + 2 * nu * a * b);
209 }
210 }
211 for (let cj = 0; cj < ny - 1; cj++) {
212 for (let ci = 0; ci < nx - 1; ci++) {
213 if (!insideCell(g, ci, cj)) continue;
214 const c = (w[id(ci + 1, cj + 1)] - w[id(ci + 1, cj)] - w[id(ci, cj + 1)] + w[id(ci, cj)]) / (hx * hy);
215 Ub += 2 * (1 - nu) * bendingStiffness(m, cellThickness(pl, ci, cj)) * c * c;
216 }
217 }
218 Ub *= 0.5 * hx * hy;
220 let Um = 0;
221 const q = (hx * hy) / 4;
222 const dNdx = new Float64Array(4), dNdy = new Float64Array(4);
223 for (let cj = 0; cj < ny - 1; cj++) {
224 for (let ci = 0; ci < nx - 1; ci++) {
225 if (!insideCell(g, ci, cj)) continue;
226 const hc = cellThickness(pl, ci, cj);
227 const C = membraneStiffness(m, hc), Gs = shearStiffness(m, hc);
228 const nn = [id(ci, cj), id(ci + 1, cj), id(ci + 1, cj + 1), id(ci, cj + 1)];
229 for (let gp = 0; gp < 4; gp++) {
230 shapeDerivatives(gp, hx, hy, dNdx, dNdy);
540490cForm the Marguerre strains from the deviation, removing a noise floor on curved designsJeremy Magland 231 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 232 for (let b = 0; b < 4; b++) {
233 const z0 = w0[nn[b]];
234 ux += dNdx[b] * u[nn[b]]; uy += dNdy[b] * u[nn[b]];
235 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 236 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 237 w0x += dNdx[b] * z0; w0y += dNdy[b] * z0;
238 }
540490cForm the Marguerre strains from the deviation, removing a noise floor on curved designsJeremy Magland 239 const exx = ux + dwx * (w0x + 0.5 * dwx);
240 const eyy = vy + dwy * (w0y + 0.5 * dwy);
241 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 242 Um += q * 0.5 * (C * (exx * exx + eyy * eyy + 2 * nu * exx * eyy) + Gs * gxy * gxy);
243 }
244 }
245 }
5778036cymbalsim: scripted plate and cymbal designer with a modal linear model and a von Karman GPU modelJeremy Magland 248}
250/**
251 * Largest angular frequency of the linearized discrete system, by power
252 * iteration on a small full grid with the same spacing and a uniform
253 * thickness h (the largest of the plate, which bounds the bending frequency;
254 * the in-plane wave speed does not depend on the thickness). The extreme
255 * mode is a grid-scale mode, so neither the plate size nor its outline
256 * matters.
257 */
258export function maxAngularFrequency(g: Grid, m: Material, h: number, nSmall = 12, iters = 300): number {
259 const gs: Grid = { nx: nSmall, ny: nSmall, hx: g.hx, hy: g.hy, n: nSmall * nSmall };
260 const pl = uniformPlate(gs, m, h);
261 const n = gs.n;
262 const x = fields(n), y = fields(n);
263 let seed = 12345;
264 const rnd = () => {
265 seed = (seed * 1664525 + 1013904223) >>> 0;
266 return seed / 4294967296 - 0.5;
267 };
268 for (let k = 0; k < n; k++) {
269 x.u[k] = rnd(); x.v[k] = rnd(); x.w[k] = rnd();
270 }
271 const eps = 1e-9;
272 let lambda = 0;
273 for (let it = 0; it < iters; it++) {
274 let norm = 0;
275 for (let k = 0; k < n; k++) norm += x.u[k] ** 2 + x.v[k] ** 2 + x.w[k] ** 2;
276 norm = Math.sqrt(norm);
277 for (let k = 0; k < n; k++) {
278 x.u[k] *= eps / norm; x.v[k] *= eps / norm; x.w[k] *= eps / norm;
279 }
5778036cymbalsim: scripted plate and cymbal designer with a modal linear model and a von Karman GPU modelJeremy Magland 281 let dot = 0;
282 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];
283 lambda = -dot / (eps * eps);
284 for (let k = 0; k < n; k++) {
285 x.u[k] = -y.u[k]; x.v[k] = -y.v[k]; x.w[k] = -y.w[k];
286 }
287 }
288 return Math.sqrt(Math.max(lambda, 0));
289}