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