1// In-browser check of the WGSL step kernel against the float64 reference:
2// one step from a smooth random state of a few thicknesses amplitude, with
3// the losses and the stand (lift-off and the felt's chain included),
4// compared as accelerations, and the update of the felt's state.
6import { type Plate, type SimParams, meanThickness, nodeX, nodeY, insideNode, standRates } from "../params.ts";
7import { accelerations, addDamping, addPadMemory, fields, fullLoss, padState, stepPadMemory } from "../model/reference.ts";
8import { fitSigma, inPlaneRate } from "../loss.ts";
9import { Simulation } from "./simulation.ts";
11export interface SelfTestResult {
12 n: number;
13 /** Largest acceleration error over all nodes and components relative to the largest reference acceleration. */
14 maxRel: number;
15 /** Largest error of the felt's arm extensions after the step relative to their largest value (0 on a free plate). */
16 padRel: number;
17 dt: number;
18}
20/** A smooth random field of the given amplitude on the plate nodes (zero outside). */
21export function smoothField(pl: Plate, amplitude: number, seed: number): Float64Array {
22 const g = pl.grid;
23 const f = new Float64Array(g.n);
24 let s = seed >>> 0;
25 const rnd = () => { s = (s * 1664525 + 1013904223) >>> 0; return s / 4294967296; };
26 const terms = Array.from({ length: 6 }, () => ({ kx: (rnd() - 0.5) * 60, ky: (rnd() - 0.5) * 60, ph: rnd() * 2 * Math.PI, a: rnd() }));
27 for (let j = 0; j < g.ny; j++) {
28 for (let i = 0; i < g.nx; i++) {
29 if (!insideNode(g, i, j)) continue;
30 const x = nodeX(g, i), y = nodeY(g, j);
31 let v = 0;
32 for (const t of terms) v += t.a * Math.sin(t.kx * x + t.ky * y + t.ph);
33 f[j * g.nx + i] = v * amplitude / 3;
34 }
35 }
36 return f;
37}
39export async function runSelfTest(device: GPUDevice, pl: Plate, params: SimParams): Promise<SelfTestResult> {
40 const g = pl.grid, n = g.n;
41 const hm = meanThickness(pl);
42 const fit = fitSigma(pl);
43 const sr = standRates(pl);
44 const loss = fullLoss(fit, inPlaneRate(pl));
45 const sim = new Simulation(device, pl, params, loss);
46 const amp = 2.5 * hm;
47 const f = { u: smoothField(pl, 0.02 * amp, 1), v: smoothField(pl, 0.02 * amp, 2), w: smoothField(pl, amp, 3) };
48 const vel = { u: smoothField(pl, 0.05, 4), v: smoothField(pl, 0.05, 5), w: smoothField(pl, 1.0, 6) };
49 const pos32 = new Float32Array(4 * n), vel32 = new Float32Array(3 * n);
50 for (let k = 0; k < n; k++) {
51 pos32[k] = f.u[k]; pos32[n + k] = f.v[k]; pos32[2 * n + k] = f.w[k];
52 vel32[k] = vel.u[k]; vel32[n + k] = vel.v[k]; vel32[2 * n + k] = vel.w[k];
53 }
54 // The felt's arms start with smooth random extensions of a fraction of the amplitude.
55 const pad = padState(pl, sr);
56 const nArms = sr ? sr.chain.tau.length : 0;
57 for (let a = 0; a < 3 * nArms; a++) pad.set(smoothField(pl, 0.3 * amp, 20 + a), a * n);
58 const pad32 = new Float32Array(pad);
59 // Use the float32-rounded state as the reference input.
60 for (let k = 0; k < n; k++) {
61 f.u[k] = pos32[k]; f.v[k] = pos32[n + k]; f.w[k] = pos32[2 * n + k];
62 vel.u[k] = vel32[k]; vel.v[k] = vel32[n + k]; vel.w[k] = vel32[2 * n + k];
63 }
64 pad.set(pad32);
65 sim.writeState(pos32, vel32, pad32);
66 await sim.stepOnce();
67 const st = await sim.readState();
68 const ref = fields(n);
69 accelerations(pl, f, ref, params.padLiftoff, sr);
70 addDamping(pl, loss, vel, ref);
71 addPadMemory(pl, sr, f, pad, params.padLiftoff, ref);
72 let maxRef = 0, maxErr = 0;
73 const dt = sim.dt;
74 for (let k = 0; k < n; k++) {
75 const au = (st.vel[k] - vel.u[k]) / dt, av = (st.vel[n + k] - vel.v[k]) / dt, aw = (st.vel[2 * n + k] - vel.w[k]) / dt;
76 maxRef = Math.max(maxRef, Math.abs(ref.u[k]), Math.abs(ref.v[k]), Math.abs(ref.w[k]));
77 maxErr = Math.max(maxErr, Math.abs(au - ref.u[k]), Math.abs(av - ref.v[k]), Math.abs(aw - ref.w[k]));
78 }
79 // The felt's state after the step, from the GPU's own new velocities.
80 let padRef = 0, padErr = 0;
81 if (sr) {
82 const dw = fields(n);
83 for (let k = 0; k < n; k++) { dw.u[k] = dt * st.vel[k]; dw.v[k] = dt * st.vel[n + k]; dw.w[k] = dt * st.vel[2 * n + k]; }
84 stepPadMemory(pl, sr, dt, dw, f, pad, params.padLiftoff);
85 for (let i = 0; i < pad.length; i++) { padRef = Math.max(padRef, Math.abs(pad[i])); padErr = Math.max(padErr, Math.abs(pad[i] - st.pad[i])); }
86 }
87 sim.destroy();
88 return { n, maxRel: maxErr / maxRef, padRel: padRef > 0 ? padErr / padRef : 0, dt };
89}