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 included), compared as accelerations.
5import { type Plate, type SimParams, meanThickness, nodeX, nodeY, insideNode, standRates } from "../params.ts";
6import { accelerations, addDamping, fields, fullLoss } from "../model/reference.ts";
7import { fitSigma, decayRate } from "../loss.ts";
8import { Simulation } from "./simulation.ts";
10export interface SelfTestResult {
11 n: number;
12 /** Largest acceleration error over all nodes and components relative to the largest reference acceleration. */
13 maxRel: number;
14 dt: number;
15}
17/** A smooth random field of the given amplitude on the plate nodes (zero outside). */
18export function smoothField(pl: Plate, amplitude: number, seed: number): Float64Array {
19 const g = pl.grid;
20 const f = new Float64Array(g.n);
21 let s = seed >>> 0;
22 const rnd = () => { s = (s * 1664525 + 1013904223) >>> 0; return s / 4294967296; };
23 const terms = Array.from({ length: 6 }, () => ({ kx: (rnd() - 0.5) * 60, ky: (rnd() - 0.5) * 60, ph: rnd() * 2 * Math.PI, a: rnd() }));
24 for (let j = 0; j < g.ny; j++) {
25 for (let i = 0; i < g.nx; i++) {
26 if (!insideNode(g, i, j)) continue;
27 const x = nodeX(g, i), y = nodeY(g, j);
28 let v = 0;
29 for (const t of terms) v += t.a * Math.sin(t.kx * x + t.ky * y + t.ph);
30 f[j * g.nx + i] = v * amplitude / 3;
31 }
32 }
33 return f;
34}
36export async function runSelfTest(device: GPUDevice, pl: Plate, params: SimParams): Promise<SelfTestResult> {
37 const g = pl.grid, n = g.n;
38 const hm = meanThickness(pl);
39 const fit = fitSigma(pl.loss, pl.mat, hm);
40 const sr = standRates(pl);
41 const loss = fullLoss(fit, sr, decayRate(pl.loss, 5000));
42 const sim = new Simulation(device, pl, params, loss);
43 const amp = 2.5 * hm;
44 const f = { u: smoothField(pl, 0.02 * amp, 1), v: smoothField(pl, 0.02 * amp, 2), w: smoothField(pl, amp, 3) };
45 const vel = { u: smoothField(pl, 0.05, 4), v: smoothField(pl, 0.05, 5), w: smoothField(pl, 1.0, 6) };
46 const pos32 = new Float32Array(4 * n), vel32 = new Float32Array(3 * n);
47 for (let k = 0; k < n; k++) {
48 pos32[k] = f.u[k]; pos32[n + k] = f.v[k]; pos32[2 * n + k] = f.w[k];
49 vel32[k] = vel.u[k]; vel32[n + k] = vel.v[k]; vel32[2 * n + k] = vel.w[k];
50 }
51 // Use the float32-rounded state as the reference input.
52 for (let k = 0; k < n; k++) {
53 f.u[k] = pos32[k]; f.v[k] = pos32[n + k]; f.w[k] = pos32[2 * n + k];
54 vel.u[k] = vel32[k]; vel.v[k] = vel32[n + k]; vel.w[k] = vel32[2 * n + k];
55 }
56 sim.writeState(pos32, vel32);
57 await sim.stepOnce();
58 const st = await sim.readState();
59 const ref = fields(n);
60 accelerations(pl, f, ref, params.padLiftoff, sr);
61 addDamping(pl, loss, vel, ref, sr);
62 let maxRef = 0, maxErr = 0;
63 const dt = sim.dt;
64 for (let k = 0; k < n; k++) {
65 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;
66 maxRef = Math.max(maxRef, Math.abs(ref.u[k]), Math.abs(ref.v[k]), Math.abs(ref.w[k]));
67 maxErr = Math.max(maxErr, Math.abs(au - ref.u[k]), Math.abs(av - ref.v[k]), Math.abs(aw - ref.w[k]));
68 }
69 sim.destroy();
70 return { n, maxRel: maxErr / maxRef, dt };
71}