1// Reading and writing .cymsim files (HDF5) with h5wasm. The layout is specified in
2// cymsim/docs/format.md; this module writes exactly what the Python package writes, apart from
3// the `state` group, which holds this application's own integrator state (writer "cymsimweb").
4// The h5wasm module is passed in by the caller: the browser uses the in-memory build ("h5wasm")
5// and writes the file image to a Blob, the tests use "h5wasm/node" on the real file system.
6import type { SimulationConfig } from "../config/schema.ts";
7import type { BackendState } from "../model/simulator.ts";
8import { AccelerationStore } from "../analysis/acceleration.ts";
10export const FORMAT_VERSION = 1;
11export const STATE_WRITER = "cymsimweb";
12const ENERGY_COLUMNS = "time kinetic membrane bending exciters max_w_over_h";
13const HIT_FIELDS = ["time", "peak_force", "contact_time", "rebound_velocity", "energy_in", "energy_out"] as const;
15// h5wasm's types are loose; we use the few calls we need through this minimal shape
16export interface H5Module { File: new (path: string, mode: string) => any; ready: Promise<{ FS: any }> }
18export interface HitRecord { time: number; kind: string; peak_force: number; contact_time: number; rebound_velocity: number; energy_in: number; energy_out: number }
20/** Everything a .cymsim file holds, in memory. */
21export interface CymsimRun {
22 yaml: string; config: SimulationConfig;
23 mesh: { xy: Float64Array; z0: Float64Array; thickness: Float64Array; normal: Float64Array; triangles: Int32Array | Uint32Array; area: Float64Array; boundary: Int32Array | Uint32Array; nNodes: number; nElements: number };
24 listener: { rFront: Float64Array; rBack: Float64Array };
25 fs: number; dt: number; substeps: number; duration: number; nSamples: number; nStepsDone: number; complete: boolean; wallTime: number; decimationSpan: number;
26 audioBaffled: Float32Array; audioUnbaffled: Float32Array; // the two listener signals, final samples so far
27 acceleration: AccelerationStore | null; taps: Float64Array | null; // taps: (W, m)
28 displacement: { rate: number; times: number[]; frames: Float32Array[] } | null;
29 energy: number[][]; // rows of six
30 hits: HitRecord[];
31 state: BackendState | null; // this application's state, when written by it
32 creator: string; created: string;
33}
35export const dtypeOf = (a: ArrayLike<number>) => a instanceof Float64Array ? "<d" : a instanceof Float32Array ? "<f" : a instanceof Int16Array ? "<h" : a instanceof Uint32Array ? "<I" : "<i";
37function attr(g: any, name: string, value: number | string | Float64Array, dtype?: string) {
38 if (typeof value === "string") g.create_attribute(name, value);
39 else if (typeof value === "number") g.create_attribute(name, value, null, dtype ?? (Number.isInteger(value) ? "<i" : "<d"));
40 else g.create_attribute(name, value, [value.length], "<d");
41}
43/** Create a resizable time-major dataset (rows, cols) and fill it from blocks without holding a copy of everything. */
44function writeSeries(g: any, name: string, dtype: string, cols: number, rows: number, chunkRows: number, fill: (write: (r0: number, block: ArrayLike<number>) => void) => void) {
45 const shape = cols > 0 ? [rows, cols] : [rows], chunks = cols > 0 ? [Math.max(1, Math.min(chunkRows, rows || chunkRows)), cols] : [Math.max(1, Math.min(chunkRows, rows || chunkRows))];
46 const empty = dtype === "<h" ? new Int16Array(0) : dtype === "<d" ? new Float64Array(0) : new Float32Array(0);
47 const ds = g.create_dataset({ name, data: empty, shape: cols > 0 ? [0, cols] : [0], dtype, maxshape: cols > 0 ? [null, cols] : [null], chunks });
48 if (rows > 0) {
49 ds.resize(shape);
50 fill((r0, block) => {
51 const n = cols > 0 ? block.length / cols : block.length;
52 ds.write_slice(cols > 0 ? [[r0, r0 + n], [0, cols]] : [[r0, r0 + n]], block);
53 });
54 }
55 return ds;
56}
58/** Write ``run`` to ``path`` in h5wasm's file system (in memory in the browser, the real disk under Node). */
59export function writeCymsim(h5: H5Module, path: string, run: CymsimRun, creator: string): void {
60 const f = new h5.File(path, "w");
61 try {
62 const N = run.mesh.nNodes, E = run.mesh.nElements;
63 attr(f, "format", "cymsim"); attr(f, "format_version", FORMAT_VERSION); attr(f, "creator", creator);
64 attr(f, "created", new Date().toISOString().replace(/\.\d{3}Z$/, "+00:00"));
65 f.create_group("config").create_dataset({ name: "yaml", data: run.yaml });
67 const mesh = f.create_group("mesh");
68 attr(mesh, "n_nodes", N); attr(mesh, "n_elements", E);
69 mesh.create_dataset({ name: "xy", data: run.mesh.xy, shape: [N, 2], dtype: "<d" });
70 mesh.create_dataset({ name: "z0", data: run.mesh.z0, shape: [N], dtype: "<d" });
71 mesh.create_dataset({ name: "thickness", data: run.mesh.thickness, shape: [N], dtype: "<d" });
72 mesh.create_dataset({ name: "normal", data: run.mesh.normal, shape: [N, 3], dtype: "<d" });
73 mesh.create_dataset({ name: "triangles", data: Int32Array.from(run.mesh.triangles), shape: [E, 3], dtype: "<i" });
74 mesh.create_dataset({ name: "area", data: run.mesh.area, shape: [E], dtype: "<d" });
75 mesh.create_dataset({ name: "boundary", data: Int32Array.from(run.mesh.boundary), shape: [run.mesh.boundary.length], dtype: "<i" });
77 const lis = f.create_group("listener"), L = run.config.listener;
78 attr(lis, "position", Float64Array.from(L.position)); attr(lis, "c_air", L.c_air, "<d"); attr(lis, "rho_air", L.rho_air, "<d");
79 attr(lis, "back_diffraction_gain", L.back_diffraction_gain, "<d");
80 lis.create_dataset({ name: "r_front", data: run.listener.rFront, shape: [E], dtype: "<d" });
81 lis.create_dataset({ name: "r_back", data: run.listener.rBack, shape: [E], dtype: "<d" });
83 const rg = f.create_group("run");
84 attr(rg, "fs", run.fs); attr(rg, "dt", run.dt, "<d"); attr(rg, "substeps", run.substeps); attr(rg, "duration", run.duration, "<d");
85 attr(rg, "n_samples", run.nSamples); attr(rg, "decimation_span", run.decimationSpan); attr(rg, "n_steps_done", run.nStepsDone);
86 attr(rg, "complete", run.complete ? 1 : 0); attr(rg, "wall_time", run.wallTime, "<d");
88 const audio = f.create_group("audio");
89 attr(audio, "fs", run.fs); attr(audio, "units", "Pa");
90 for (const [name, data] of [["baffled", run.audioBaffled], ["unbaffled", run.audioUnbaffled]] as const) {
91 writeSeries(audio, name, "<f", 0, data.length, 65536, (write) => write(0, data));
92 }
94 if (run.acceleration && run.taps) {
95 const a = run.acceleration, g = f.create_group("acceleration");
96 attr(g, "fs", run.fs); attr(g, "component", "normal"); attr(g, "units", "m/s^2"); attr(g, "encoding", "int16 * scale[t]");
97 attr(g, "description", "acceleration of every node along the rest-surface normal (front, +z side, positive), low-pass filtered by the kernel in 'taps' and sampled at fs; sample j is centred on t = j / fs");
98 const rows = Math.max(1, Math.min(1024, Math.floor(16e6 / (2 * N))));
99 writeSeries(g, "normal", "<h", N, a.nSamples, rows, (write) => { for (const b of a.blocks) write(b.first, b.data); });
100 writeSeries(g, "scale", "<f", 0, a.nSamples, 65536, (write) => { for (const b of a.blocks) write(b.first, b.scale); });
101 const W = run.decimationSpan, m = run.substeps;
102 const taps = g.create_dataset({ name: "taps", data: run.taps, shape: [W, m], dtype: "<d" });
103 attr(taps, "layout", "taps[w, kk]: weight of sub-step j0 m + kk in audio sample j0 + w - W/2 + 1");
104 attr(taps, "kaiser_beta", 8.0, "<d"); attr(taps, "cutoff", 0.46, "<d");
105 }
107 if (run.displacement) {
108 const d = run.displacement, g = f.create_group("displacement");
109 attr(g, "rate", d.rate, "<d"); attr(g, "units", "m"); attr(g, "description", "transverse displacement w of every node at the frame times");
110 writeSeries(g, "w", "<f", N, d.frames.length, 64, (write) => d.frames.forEach((fr, i) => write(i, fr)));
111 writeSeries(g, "time", "<d", 0, d.times.length, 4096, (write) => write(0, Float64Array.from(d.times)));
112 }
114 const eg = f.create_group("energy");
115 const table = writeSeries(eg, "table", "<d", 6, run.energy.length, 1024, (write) => write(0, Float64Array.from(run.energy.flat())));
116 attr(table, "columns", ENERGY_COLUMNS);
118 const hg = f.create_group("hits"), K = run.hits.length;
119 for (const name of HIT_FIELDS) hg.create_dataset({ name, data: Float64Array.from(run.hits, (h) => h[name]), shape: [K], dtype: "<d" });
120 hg.create_dataset({ name: "kind", data: run.hits.map((h) => h.kind), shape: [K], dtype: "S16" });
122 if (run.state) {
123 const s = run.state, g = f.create_group("state");
124 attr(g, "writer", STATE_WRITER); attr(g, "n_steps_done", s.step); attr(g, "samples_done", s.samplesDone); attr(g, "backend", s.kind);
125 for (const [name, a] of Object.entries(s.arrays)) g.create_dataset({ name, data: a, shape: [a.length], dtype: dtypeOf(a) });
126 }
127 } finally { f.close(); }
128}
130const str = (v: unknown) => (typeof v === "string" ? v : Array.isArray(v) ? String(v[0]) : String(v));
131const num = (v: unknown) => Number(Array.isArray(v) ? v[0] : v);
133/** Read a .cymsim file from h5wasm's file system. ``parse`` turns the YAML text into a configuration. */
134export function readCymsim(h5: H5Module, path: string, parse: (yaml: string) => SimulationConfig): CymsimRun {
135 const f = new h5.File(path, "r");
136 try {
137 if (str(f.get_attribute("format", true)) !== "cymsim") throw new Error("not a .cymsim file");
138 const get = (p: string) => { const d = f.get(p); if (!d) throw new Error(`missing ${p}`); return d; };
139 const has = (p: string) => f.get(p) !== null && f.get(p) !== undefined;
140 const value = <T,>(p: string) => get(p).value as T;
141 const a = (p: string, name: string) => get(p).get_attribute(name, true);
142 const yamlText = str(value("config/yaml"));
143 const config = parse(yamlText);
144 const nNodes = num(a("mesh", "n_nodes")), nElements = num(a("mesh", "n_elements"));
145 const f64 = (p: string) => Float64Array.from(value<ArrayLike<number>>(p));
146 const f32 = (p: string) => Float32Array.from(value<ArrayLike<number>>(p));
147 const fs = num(a("run", "fs")), substeps = num(a("run", "substeps"));
148 let acceleration: AccelerationStore | null = null, taps: Float64Array | null = null;
149 if (has("acceleration")) {
150 acceleration = new AccelerationStore(nNodes, fs);
151 const ds = get("acceleration/normal"), scale = f32("acceleration/scale"), T = ds.shape[0] as number;
152 // read in slabs of about 64 MB so that no second copy of the whole record is needed
153 const slab = Math.max(1, Math.floor(32e6 / nNodes));
154 for (let r0 = 0; r0 < T; r0 += slab) {
155 const r1 = Math.min(T, r0 + slab);
156 const data = ds.slice([[r0, r1], [0, nNodes]]) as Int16Array;
157 acceleration.pushEncoded(Int16Array.from(data), scale.slice(r0, r1));
158 }
159 taps = f64("acceleration/taps");
160 }
161 let displacement: CymsimRun["displacement"] = null;
162 if (has("displacement")) {
163 const w = get("displacement/w"), Td = w.shape[0] as number, all = w.value as Float32Array, times = Array.from(f64("displacement/time"));
164 const frames: Float32Array[] = [];
165 for (let i = 0; i < Td; i++) frames.push(Float32Array.from(all.subarray(i * nNodes, (i + 1) * nNodes)));
166 displacement = { rate: num(a("displacement", "rate")), times, frames };
167 }
168 const et = get("energy/table"), erows = et.shape[0] as number, eflat = Array.from(et.value as ArrayLike<number>);
169 const energy: number[][] = [];
170 for (let i = 0; i < erows; i++) energy.push(eflat.slice(6 * i, 6 * i + 6));
171 const K = get("hits/time").shape[0] as number, kinds = get("hits/kind").value as unknown as string[];
172 const hits: HitRecord[] = [];
173 for (let k = 0; k < K; k++) {
174 const h: any = { kind: str(Array.isArray(kinds) ? kinds[k] : kinds) };
175 for (const name of HIT_FIELDS) h[name] = Number((get(`hits/${name}`).value as ArrayLike<number>)[k]);
176 hits.push(h);
177 }
178 let state: BackendState | null = null;
179 if (has("state") && str(a("state", "writer")) === STATE_WRITER) {
180 const g = get("state"), arrays: BackendState["arrays"] = {};
181 for (const name of g.keys() as string[]) {
182 const v = g.get(name).value;
183 arrays[name] = v instanceof Float64Array || v instanceof Float32Array || v instanceof Uint32Array ? v : Float32Array.from(v as ArrayLike<number>);
184 }
185 state = { kind: str(a("state", "backend")) as "cpu" | "gpu", step: num(a("state", "n_steps_done")), samplesDone: num(a("state", "samples_done")), arrays };
186 }
187 const optional = (name: string, dflt: number) => { try { return num(a("run", name)); } catch { return dflt; } };
188 return {
189 yaml: yamlText, config,
190 mesh: { xy: f64("mesh/xy"), z0: f64("mesh/z0"), thickness: f64("mesh/thickness"), normal: f64("mesh/normal"), triangles: Int32Array.from(value<ArrayLike<number>>("mesh/triangles")),
191 area: f64("mesh/area"), boundary: Int32Array.from(value<ArrayLike<number>>("mesh/boundary")), nNodes, nElements },
192 listener: { rFront: f64("listener/r_front"), rBack: f64("listener/r_back") },
193 fs, dt: num(a("run", "dt")), substeps, duration: num(a("run", "duration")), nSamples: num(a("run", "n_samples")),
194 nStepsDone: optional("n_steps_done", 0), complete: optional("complete", 0) !== 0, wallTime: optional("wall_time", 0), decimationSpan: optional("decimation_span", 32),
195 audioBaffled: f32("audio/baffled"), audioUnbaffled: f32("audio/unbaffled"), acceleration, taps, displacement, energy, hits, state,
196 creator: str(f.get_attribute("creator", true)), created: str(f.get_attribute("created", true)),
197 };
198 } finally { f.close(); }
199}