Steps the fluid's velocity and density on a grid of cells to a steady state and sums the mass through openings and the force on walls.
How the precision is kept
The populations, densities and velocities are Float64Array and each step is 64-bit arithmetic. The result's accuracy is set by the method: the cell size, the lattice speed chosen to keep density changes to a few percent, the wall treatment, and the steady-state tolerance. Fluid properties are table values rounded to 4 or 5 digits.
Calculations under this policy
10 calculations in 1 source file. The code is read from the released app when this page loads. It is the calculation people run, not a copy kept for this page.
Air density is computed from 101325 ÷ (287.05 × 293.15) and then rounded to 4 decimals; air viscosity from Sutherland's law and then rounded to 4 significant digits. Water uses the two table values 998.21 and 1.0016e-3. Every flow result is proportional to these rounded inputs.
Connects each fluid cell to its neighbors and sets the term of every boundary link.
moving wall: a = 6 wᵢ (cᵢ · u); interpolated wall when |q − ½| > 1e-6, q clamped to [0.001, 1]
Link terms and wall distances are stored in Float64Array. The wall-surface normals kept for the turbulence model are stored with Float32Array.from; the laminar solver does not read them.
flow.js
lines 56 to 136
PartCraft 0.8.1
/**
* A lattice from a cell classification.
* dom = { nx, ny, nz, kind: Int32Array(nx·ny·nz): −1 solid/outside, 0 fluid, g ≥ 1 a boundary cell of group g,
* periodic: [bx, by, bz], groups: [null, {type 'wall'|'velocity'|'pressure'|'slip', u?:[3], rho?, axis?, force?}] (index ≥ 1),
* wallGroup(cell, i) → {g, q, d?, n?} | null : what a link from fluid cell `cell` along direction i (towards a solid cell)
* hits: a group (≥ 1, e.g. an inlet drawn on a capped face) and the wall distance q in (0, 1] used by the link (null =
* plain wall, q ½); d, n: the true distance fraction and surface normal (into the solid) for the turbulence wall model,
* profile(g, x, y, z) → [ux, uy, uz] | null : a velocity group's value at a boundary cell (tests: parabolic inflow) }
* Wall group 0 is the part (counted for forces). Returns the lattice (arrays) — solve with `run`.
*/
export function buildLattice(dom) {
const { nx, ny, nz, kind } = dom, per = dom.periodic || [false, false, false], groups = dom.groups || [null];
const ncell = nx * ny * nz;
const index = new Int32Array(ncell).fill(-1);
let nf = 0;
for (let c = 0; c < ncell; c++) if (kind[c] === 0) index[c] = nf++;
if (!nf) throw new Error(t('study.flow.err.no_fluid_space'));
const cells = new Int32Array(nf);
for (let c = 0; c < ncell; c++) if (index[c] >= 0) cells[index[c]] = c;
const nbo = new Int32Array(nf * Q);
const bulk = new Uint8Array(nf);
const L = { kind: [], cell: [], dir: [], a: [], src: [], g: [], force: [], wd: [], wn: [] };
const wrap = (v, n, p) => (v >= 0 && v < n ? v : p ? (v + n) % n : -1);
const cellAt = (i, j, k) => (i < 0 || j < 0 || k < 0 ? -1 : i + nx * (j + ny * k));
for (let n = 0; n < nf; n++) {
const c = cells[n], i0 = c % nx, j0 = ((c - i0) / nx) % ny, k0 = (c - i0 - nx * j0) / (nx * ny);
let isBulk = 1;
nbo[n * Q] = n * Q;
for (let q = 1; q < Q; q++) {
// pull: the population moving along c_q arrives from x − c_q
const si = wrap(i0 - CX[q], nx, per[0]), sj = wrap(j0 - CY[q], ny, per[1]), sk = wrap(k0 - CZ[q], nz, per[2]);
const sc = cellAt(si, sj, sk), sn = sc >= 0 ? index[sc] : -1;
if (sn >= 0) { nbo[n * Q + q] = sn * Q + q; continue; }
isBulk = 0;
const sg = sc >= 0 ? kind[sc] : -1, oq = OPP[q];
let lk = K_WALL, a = 0.5, src = -1, g = 0, force = 1, wd = NaN, wn = null;
if (sg >= 1) {
g = sg; const G = groups[g];
if (!G) throw new Error('A boundary cell names no group.'); // i18n-ignore: internal check, never reaches a person
if (G.type === 'velocity') {
const u = (dom.profile && dom.profile(g, i0 - CX[q], j0 - CY[q], k0 - CZ[q])) || G.u;
lk = K_VEL; a = 6 * W[q] * (CX[q] * u[0] + CY[q] * u[1] + CZ[q] * u[2]); force = 0;
} else if (G.type === 'pressure') { lk = K_PRES; a = G.rho; force = 0; }
else if (G.type === 'slip') {
lk = K_SLIP; force = 0;
// specular: the population arriving along c_q left x − (tangential part of c_q) along c_q mirrored in the normal axis
const ax = G.axis, t = [CX[q], CY[q], CZ[q]], m = t.slice(); t[ax] = 0; m[ax] = -m[ax];
const mi = findDir(m[0], m[1], m[2]);
const xi = wrap(i0 - t[0], nx, per[0]), xj = wrap(j0 - t[1], ny, per[1]), xk = wrap(k0 - t[2], nz, per[2]);
const xc = cellAt(xi, xj, xk), xn = xc >= 0 ? index[xc] : -1;
src = xn >= 0 && mi >= 0 ? xn * Q + mi : -1;
} else { lk = K_WALL; force = G.force ? 1 : 0; }
} else {
// a solid (or outside) cell: the outgoing population ō = OPP[q] travels x → x + c_ō = x − c_q and hits a wall there
const wg = sc >= 0 && dom.wallGroup ? dom.wallGroup(c, oq) : null;
if (sc < 0) force = 0; // off the grid (tests' channel walls): a wall, not the part
if (wg) { wd = wg.d ?? wg.q ?? NaN; wn = wg.n || null; }
if (wg && wg.g >= 1) {
g = wg.g; const G = groups[g];
if (G.type === 'velocity') { lk = K_VEL; a = 6 * W[q] * (CX[q] * G.u[0] + CY[q] * G.u[1] + CZ[q] * G.u[2]); force = 0; }
else if (G.type === 'pressure') { lk = K_PRES; a = G.rho; force = 0; }
else { lk = K_WALL; force = G.force ? 1 : 0; }
} else if (wg && wg.q !== undefined && Math.abs(wg.q - 0.5) > 1e-6) {
lk = K_BOUZIDI; a = Math.min(1, Math.max(1e-3, wg.q));
// the second fluid node behind x (x + c_q) for q < ½
const fi = wrap(i0 + CX[q], nx, per[0]), fj = wrap(j0 + CY[q], ny, per[1]), fk = wrap(k0 + CZ[q], nz, per[2]);
const fc = cellAt(fi, fj, fk), fn = fc >= 0 ? index[fc] : -1;
src = fn >= 0 ? fn * Q + oq : -1;
if (a < 0.5 && src < 0) { lk = K_WALL; a = 0.5; }
}
}
nbo[n * Q + q] = nf * Q + L.kind.length;
L.kind.push(lk); L.cell.push(n); L.dir.push(q); L.a.push(a); L.src.push(src); L.g.push(g); L.force.push(force);
L.wd.push(wd); L.wn.push(...(wn || [NaN, NaN, NaN]));
}
bulk[n] = isBulk;
}
const links = { kind: Uint8Array.from(L.kind), cell: Int32Array.from(L.cell), dir: Uint8Array.from(L.dir), a: Float64Array.from(L.a), src: Int32Array.from(L.src), g: Int32Array.from(L.g), force: Uint8Array.from(L.force), n: L.kind.length,
wd: Float64Array.from(L.wd), wn: Float32Array.from(L.wn) }; // wall links: true distance fraction and surface normal (NaN: not known)
return { nx, ny, nz, nf, cells, index, nbo, bulk, links, groups, ncell, periodic: per.slice() };
}
64-bit. The function refuses a relaxation time that is not above ½, so the division by τ − ½ cannot be by zero; a value just above ½ gives a very large τ⁻.
flow.js
lines 160 to 178
PartCraft 0.8.1
/**
* A solver state on a lattice: tau (symmetric relaxation time, ν = (τ − ½)/3), initial velocity u0 (lattice).
* Options for turbulence (flow-turb.js) and tests: eddy: true gives the state a per-cell eddy viscosity `nut` (lattice,
* Float64Array(nf), 0 at first) and the eddy collision (collideEddy); with it every wall link joins the local mass
* correction (wall-model slip terms must pass no mass). smagorinsky: the stabilizer's constant (eddy only, 0 = off).
* force: [fx, fy, fz] a uniform body force per unit volume (lattice; added as (1 − ω−/2)·3 w_i c_i·F after collision,
* velocity u = (Σ c f + F/2)/ρ) — periodic test channels.
*/
export function makeState(lat, { tau, u0 = [0, 0, 0], magic = MAGIC, eddy = false, force = null, smagorinsky = 0 }) {
if (!(tau > 0.5)) throw new Error('The relaxation time must be above one half.'); // i18n-ignore: caller error (scaling() refuses first)
const nf = lat.nf, A = new Float64Array(nf * Q + lat.links.n), B = new Float64Array(nf * Q + lat.links.n);
const ux = new Float64Array(nf).fill(u0[0]), uy = new Float64Array(nf).fill(u0[1]), uz = new Float64Array(nf).fill(u0[2]), rho = new Float64Array(nf).fill(1);
const usq = u0[0] ** 2 + u0[1] ** 2 + u0[2] ** 2;
for (let n = 0; n < nf; n++) for (let q = 0; q < Q; q++) { const cu = CX[q] * u0[0] + CY[q] * u0[1] + CZ[q] * u0[2]; A[n * Q + q] = W[q] * (1 + 3 * cu + 4.5 * cu * cu - 1.5 * usq); }
const tm = 0.5 + magic / (tau - 0.5);
const S = { lat, src: A, dst: B, ux, uy, uz, rho, tau, tauMinus: tm, wp: 1 / tau, wm: 1 / tm, magic, step: 0, smagorinsky, omegaMinus: OMEGA_MINUS, nueff: null, wallSlot: null, wallStress: null, nut: eddy ? new Float64Array(nf) : null, force: force && force.some((v) => v) ? Float64Array.from(force) : null, turb: null, ...linearLinks(lat, eddy) };
S.whole = slab(S, 0, nf);
return S;
}
Boundary links as linear terms
64-bit
Writes every wall, moving-wall and slip link as two coefficients and a constant, for the time step's inner loop.
The coefficients are computed once in 64-bit and stored in Float64Array. They are the same expressions as in the per-link function incoming. In turbulent runs the wall model later rewrites the constant a of wall links.
flow.js
lines 180 to 201
PartCraft 0.8.1
/** Every boundary link but the pressure ones as f = c1·src[o1] + c2·src[o2] + a (bounce-back, Bouzidi, slip, moving wall), for a tight loop. */
function linearLinks(lat, allWalls = false) {
const L = lat.links, k = [], o1 = [], o2 = [], c1 = [], c2 = [], a = [], prs = [];
for (let i = 0; i < L.n; i++) {
const b = L.cell[i] * Q, q = L.dir[i], own = b + OPP[q];
switch (L.kind[i]) {
case K_PRES: prs.push(i); continue;
case K_BOUZIDI: {
const qq = L.a[i];
if (qq < 0.5) { k.push(i); o1.push(own); o2.push(L.src[i]); c1.push(2 * qq); c2.push(1 - 2 * qq); a.push(0); }
else { k.push(i); o1.push(own); o2.push(b + q); c1.push(1 / (2 * qq)); c2.push((2 * qq - 1) / (2 * qq)); a.push(0); }
continue;
}
case K_VEL: k.push(i); o1.push(own); o2.push(own); c1.push(1); c2.push(0); a.push(L.a[i]); continue;
case K_SLIP: k.push(i); o1.push(L.src[i] >= 0 ? L.src[i] : own); o2.push(own); c1.push(1); c2.push(0); a.push(0); continue;
default: k.push(i); o1.push(own); o2.push(own); c1.push(1); c2.push(0); a.push(0);
}
}
const bz = { k: [], cell: [], own: [] };
if (lat.massFix !== false) for (let i = 0; i < L.n; i++) if (L.kind[i] === K_BOUZIDI || (allWalls && L.kind[i] === K_WALL)) { bz.k.push(i); bz.cell.push(L.cell[i]); bz.own.push(L.cell[i] * Q + OPP[L.dir[i]]); }
return { bz: { k: Int32Array.from(bz.k), cell: Int32Array.from(bz.cell), own: Int32Array.from(bz.own) }, corr: new Float64Array(lat.nf), lin: { k: Int32Array.from(k), cell: Int32Array.from(k, (i) => L.cell[i]), o1: Int32Array.from(o1), o2: Int32Array.from(o2), c1: Float64Array.from(c1), c2: Float64Array.from(c2), a: Float64Array.from(a) }, prs: Int32Array.from(prs), prsCell: Int32Array.from(prs, (i) => L.cell[i]) };
}
One time step: stream and collide
64-bit
Pulls the populations into each cell, takes density and velocity, and relaxes toward equilibrium.
ρ = Σ fᵢ; u = (Σ cᵢ fᵢ + F ÷ 2) ÷ ρ
Density is a 19-term sum in a fixed order; velocity is a sum of differences of opposite populations times 1 ÷ ρ. Populations are of the order of 0.03 to 0.3, so the velocity's absolute rounding error is of the order of 1e-17 in lattice units whatever the speed, and slow regions keep fewer relative digits than fast ones.
flow.js
lines 236 to 314
PartCraft 0.8.1
/**
* The kernel on one slab. The boundary links' arriving populations are computed first into ghost slots after the cells'
* (src[nf·19 + k]), which only the link's own cell reads (its `nbo` entry points there) — so the cell loop has no
* branches: pull, collide (TRT; collideEddy with an eddy viscosity), store.
*/
export function stepSlab(S, sl) {
const { lat, src, dst, ux, uy, uz, rho, wp, wm, nut, force } = S, { nbo, nf } = lat, G = nf * Q;
const P = S.lin, PK = P.k, O1 = P.o1, O2 = P.o2, C1 = P.c1, C2 = P.c2, PA = P.a, PR = S.prs;
for (let m = sl.l0; m < sl.l1; m++) src[G + PK[m]] = C1[m] * src[O1[m]] + C2[m] * src[O2[m]] + PA[m];
// interpolated walls: each link's mass defect (what left minus what came back) goes back into the node's rest population
const BZ = S.bz, corr = S.corr;
if (BZ.k.length) { corr.fill(0, sl.n0, sl.n1); for (let m = sl.b0; m < sl.b1; m++) corr[BZ.cell[m]] += src[BZ.own[m]] - src[G + BZ.k[m]]; }
for (let m = sl.p0; m < sl.p1; m++) { const k = PR[m]; src[G + k] = incoming(lat, k, src, ux, uy, uz); }
if (nut) { collideEddy(S, sl); return; }
const Fx = force ? force[0] : 0, Fy = force ? force[1] : 0, Fz = force ? force[2] : 0;
let f0, f1, f2, f3, f4, f5, f6, f7, f8, f9, f10, f11, f12, f13, f14, f15, f16, f17, f18, ep, em, sp, sm, w;
for (let n = sl.n0; n < sl.n1; n++) {
const b = n * Q;
f0 = src[b] + corr[n];
f1 = src[nbo[b + 1]];
f2 = src[nbo[b + 2]];
f3 = src[nbo[b + 3]];
f4 = src[nbo[b + 4]];
f5 = src[nbo[b + 5]];
f6 = src[nbo[b + 6]];
f7 = src[nbo[b + 7]];
f8 = src[nbo[b + 8]];
f9 = src[nbo[b + 9]];
f10 = src[nbo[b + 10]];
f11 = src[nbo[b + 11]];
f12 = src[nbo[b + 12]];
f13 = src[nbo[b + 13]];
f14 = src[nbo[b + 14]];
f15 = src[nbo[b + 15]];
f16 = src[nbo[b + 16]];
f17 = src[nbo[b + 17]];
f18 = src[nbo[b + 18]];
const r = f0 + f1 + f2 + f3 + f4 + f5 + f6 + f7 + f8 + f9 + f10 + f11 + f12 + f13 + f14 + f15 + f16 + f17 + f18;
const ir = 1 / r;
const u0 = (f1 - f2 + f7 - f8 + f9 - f10 + f11 - f12 + f13 - f14 + 0.5 * Fx) * ir;
const u1 = (f3 - f4 + f7 - f8 - f9 + f10 + f15 - f16 + f17 - f18 + 0.5 * Fy) * ir;
const u2 = (f5 - f6 + f11 - f12 - f13 + f14 + f15 - f16 - f17 + f18 + 0.5 * Fz) * ir;
rho[n] = r; ux[n] = u0; uy[n] = u1; uz[n] = u2;
const us = 1.5 * (u0 * u0 + u1 * u1 + u2 * u2);
dst[b] = f0 - wp * (f0 - r * (1 / 3) * (1 - us));
// TRT on the 9 pairs (a, a+1): symmetric part relaxed with ω+, antisymmetric with ω−
w = r / 18;
ep = w * (1 + 4.5 * (u0) * (u0) - us); em = w * 3 * (u0);
sp = wp * (0.5 * (f1 + f2) - ep); sm = wm * (0.5 * (f1 - f2) - em);
dst[b + 1] = f1 - sp - sm; dst[b + 2] = f2 - sp + sm;
ep = w * (1 + 4.5 * (u1) * (u1) - us); em = w * 3 * (u1);
sp = wp * (0.5 * (f3 + f4) - ep); sm = wm * (0.5 * (f3 - f4) - em);
dst[b + 3] = f3 - sp - sm; dst[b + 4] = f4 - sp + sm;
ep = w * (1 + 4.5 * (u2) * (u2) - us); em = w * 3 * (u2);
sp = wp * (0.5 * (f5 + f6) - ep); sm = wm * (0.5 * (f5 - f6) - em);
dst[b + 5] = f5 - sp - sm; dst[b + 6] = f6 - sp + sm;
w = r / 36;
ep = w * (1 + 4.5 * (u0 + u1) * (u0 + u1) - us); em = w * 3 * (u0 + u1);
sp = wp * (0.5 * (f7 + f8) - ep); sm = wm * (0.5 * (f7 - f8) - em);
dst[b + 7] = f7 - sp - sm; dst[b + 8] = f8 - sp + sm;
ep = w * (1 + 4.5 * (u0 - u1) * (u0 - u1) - us); em = w * 3 * (u0 - u1);
sp = wp * (0.5 * (f9 + f10) - ep); sm = wm * (0.5 * (f9 - f10) - em);
dst[b + 9] = f9 - sp - sm; dst[b + 10] = f10 - sp + sm;
ep = w * (1 + 4.5 * (u0 + u2) * (u0 + u2) - us); em = w * 3 * (u0 + u2);
sp = wp * (0.5 * (f11 + f12) - ep); sm = wm * (0.5 * (f11 - f12) - em);
dst[b + 11] = f11 - sp - sm; dst[b + 12] = f12 - sp + sm;
ep = w * (1 + 4.5 * (u0 - u2) * (u0 - u2) - us); em = w * 3 * (u0 - u2);
sp = wp * (0.5 * (f13 + f14) - ep); sm = wm * (0.5 * (f13 - f14) - em);
dst[b + 13] = f13 - sp - sm; dst[b + 14] = f14 - sp + sm;
ep = w * (1 + 4.5 * (u1 + u2) * (u1 + u2) - us); em = w * 3 * (u1 + u2);
sp = wp * (0.5 * (f15 + f16) - ep); sm = wm * (0.5 * (f15 - f16) - em);
dst[b + 15] = f15 - sp - sm; dst[b + 16] = f16 - sp + sm;
ep = w * (1 + 4.5 * (u1 - u2) * (u1 - u2) - us); em = w * 3 * (u1 - u2);
sp = wp * (0.5 * (f17 + f18) - ep); sm = wm * (0.5 * (f17 - f18) - em);
dst[b + 17] = f17 - sp - sm; dst[b + 18] = f18 - sp + sm;
}
// body force on the antisymmetric part, (1 − ω−/2)·3 w c·F (Guo et al. 2002 in TRT form; consistent with the F/2 velocity shift)
if (force) for (let n = sl.n0; n < sl.n1; n++) for (let q = 1; q < Q; q += 2) { const g = (1 - 0.5 * wm) * 3 * W[q] * (CX[q] * Fx + CY[q] * Fy + CZ[q] * Fz); dst[n * Q + q] += g; dst[n * Q + q + 1] -= g; }
}
Collision with an eddy viscosity
64-bit
The collision used in turbulent runs, with a per-cell viscosity and a stabilizer.
64-bit throughout. The non-equilibrium stress is formed as population sums minus twice the equilibrium, a subtraction of nearly equal numbers when the flow is close to equilibrium; the stabilizer and the wall stress are built on that difference.
flow.js
lines 317 to 402
PartCraft 0.8.1
/**
* The collision with a per-cell eddy viscosity (turbulent runs): TRT, τ+ = ½ + 3(ν + ν_t) per cell, τ− from the magic
* Λ = 3/16 but with ω− ≥ OMEGA_MINUS, plus a STABILIZER: τ+ = max(τ+, τ_s), τ_s the Smagorinsky relaxation time of the
* resolved strain read from the non-equilibrium stress Π^neq (Hou et al. 1996 with ν → 0: τ_s = ¼ + √(1/16 + (9/2)√2
* C²|Π^neq|/ρ), C = S.smagorinsky). In an equilibrium boundary layer the model's ν_t = κu_τy exceeds C²|S| = C²u_τ/(κy)
* wherever y > C/κ cells (0.41 cells at C 0.17), so it does not act on a wall-function profile, nor in a free stream; it
* acts where the grid sees gradients the turbulence model does not cover (sharp inlet corners, start-up transients),
* which a lattice with the molecular τ at ½ cannot survive otherwise. The bulk (trace) part of Π^neq relaxes at
* BULK_OMEGA. S.nueff (if present) receives the viscosity used, ν_eff − ν; for cells with S.wallSlot[n] = w ≥ 0,
* S.wallStress[6w…] receives the viscous stress −(1 − ω+/2)Π^neq (xx yy zz xy xz yz) for the wall model.
*/
function collideEddy(S, sl) {
const { lat, src, dst, ux, uy, uz, rho, nut, force } = S, nbo = lat.nbo, nu0 = (S.tau - 0.5) / 3, corr = S.corr, mg = S.magic;
const cs2 = 4.5 * Math.SQRT2 * S.smagorinsky ** 2, nue = S.nueff, wmMin = S.omegaMinus, wsl = S.wallSlot, wst = S.wallStress;
const Fx = force ? force[0] : 0, Fy = force ? force[1] : 0, Fz = force ? force[2] : 0;
for (let n = sl.n0; n < sl.n1; n++) {
const b = n * Q;
const f0 = src[b] + corr[n];
const f1 = src[nbo[b + 1]];
const f2 = src[nbo[b + 2]];
const f3 = src[nbo[b + 3]];
const f4 = src[nbo[b + 4]];
const f5 = src[nbo[b + 5]];
const f6 = src[nbo[b + 6]];
const f7 = src[nbo[b + 7]];
const f8 = src[nbo[b + 8]];
const f9 = src[nbo[b + 9]];
const f10 = src[nbo[b + 10]];
const f11 = src[nbo[b + 11]];
const f12 = src[nbo[b + 12]];
const f13 = src[nbo[b + 13]];
const f14 = src[nbo[b + 14]];
const f15 = src[nbo[b + 15]];
const f16 = src[nbo[b + 16]];
const f17 = src[nbo[b + 17]];
const f18 = src[nbo[b + 18]];
const r = f0 + f1 + f2 + f3 + f4 + f5 + f6 + f7 + f8 + f9 + f10 + f11 + f12 + f13 + f14 + f15 + f16 + f17 + f18, ir = 1 / r;
const u0 = (f1 - f2 + f7 - f8 + f9 - f10 + f11 - f12 + f13 - f14 + 0.5 * Fx) * ir;
const u1 = (f3 - f4 + f7 - f8 - f9 + f10 + f15 - f16 + f17 - f18 + 0.5 * Fy) * ir;
const u2 = (f5 - f6 + f11 - f12 - f13 + f14 + f15 - f16 - f17 + f18 + 0.5 * Fz) * ir;
rho[n] = r; ux[n] = u0; uy[n] = u1; uz[n] = u2;
const us = 1.5 * (u0 * u0 + u1 * u1 + u2 * u2), w1 = r / 18, w2 = r / 36;
const c0 = u0, e0 = w1 * (1 + 4.5 * c0 * c0 - us), d0 = f1 + f2 - 2 * e0;
const c1 = u1, e1 = w1 * (1 + 4.5 * c1 * c1 - us), d1 = f3 + f4 - 2 * e1;
const c2 = u2, e2 = w1 * (1 + 4.5 * c2 * c2 - us), d2 = f5 + f6 - 2 * e2;
const c3 = u0 + u1, e3 = w2 * (1 + 4.5 * c3 * c3 - us), d3 = f7 + f8 - 2 * e3;
const c4 = u0 - u1, e4 = w2 * (1 + 4.5 * c4 * c4 - us), d4 = f9 + f10 - 2 * e4;
const c5 = u0 + u2, e5 = w2 * (1 + 4.5 * c5 * c5 - us), d5 = f11 + f12 - 2 * e5;
const c6 = u0 - u2, e6 = w2 * (1 + 4.5 * c6 * c6 - us), d6 = f13 + f14 - 2 * e6;
const c7 = u1 + u2, e7 = w2 * (1 + 4.5 * c7 * c7 - us), d7 = f15 + f16 - 2 * e7;
const c8 = u1 - u2, e8 = w2 * (1 + 4.5 * c8 * c8 - us), d8 = f17 + f18 - 2 * e8;
const pxx = d0 + d3 + d4 + d5 + d6, pyy = d1 + d3 + d4 + d7 + d8, pzz = d2 + d5 + d6 + d7 + d8;
const pxy = d3 - d4, pxz = d5 - d6, pyz = d7 - d8;
let t = 0.5 + 3 * (nu0 + nut[n]);
if (cs2) { const ts = 0.25 + Math.sqrt(0.0625 + cs2 * Math.sqrt(pxx * pxx + pyy * pyy + pzz * pzz + 2 * (pxy * pxy + pxz * pxz + pyz * pyz)) * ir); if (ts > t) t = ts; }
if (nue) nue[n] = (t - 0.5) / 3 - nu0;
// ω− = max(OMEGA_MINUS, magic): at the small τ+ of turbulent runs the magic Λ would leave the odd moments all but
// unrelaxed (ω− → 0), which staircase walls and corners cannot survive
const wp = 1 / t, tm = t - 0.5, wm = Math.max(wmMin, tm / (0.5 * tm + mg));
if (wsl) { const w = wsl[n]; if (w >= 0) { const k = -(1 - 0.5 * wp), o = 6 * w; wst[o] = k * pxx; wst[o + 1] = k * pyy; wst[o + 2] = k * pzz; wst[o + 3] = k * pxy; wst[o + 4] = k * pxz; wst[o + 5] = k * pyz; } }
// bulk (trace) part of Π^neq relaxed at ω_b instead of ω+: Δf_i = −(ω_b − ω+)·(9/2) w_i (|c_i|² − 1) tr Π/3 changes
// the trace only (mass, momentum and the deviatoric stress are untouched) — damps the sound waves a lattice at τ ≈ ½
// would otherwise keep ringing between inlets and outlets
const tb = (BULK_OMEGA - wp) * 1.5 * (pxx + pyy + pzz), te = tb / 36;
dst[b] = f0 - wp * (f0 - r * (1 / 3) * (1 - us)) + tb / 3;
const kf = force ? (1 - 0.5 * wm) * 3 : 0;
{ const sp = wp * (0.5 * (f1 + f2) - e0), sm = wm * (0.5 * (f1 - f2) - w1 * 3 * c0), g = kf * (Fx) * 1 / 18;
dst[b + 1] = f1 - sp - sm + g; dst[b + 2] = f2 - sp + sm - g; }
{ const sp = wp * (0.5 * (f3 + f4) - e1), sm = wm * (0.5 * (f3 - f4) - w1 * 3 * c1), g = kf * (Fy) * 1 / 18;
dst[b + 3] = f3 - sp - sm + g; dst[b + 4] = f4 - sp + sm - g; }
{ const sp = wp * (0.5 * (f5 + f6) - e2), sm = wm * (0.5 * (f5 - f6) - w1 * 3 * c2), g = kf * (Fz) * 1 / 18;
dst[b + 5] = f5 - sp - sm + g; dst[b + 6] = f6 - sp + sm - g; }
{ const sp = wp * (0.5 * (f7 + f8) - e3), sm = wm * (0.5 * (f7 - f8) - w2 * 3 * c3), g = kf * (Fx + Fy) * 1 / 36;
dst[b + 7] = f7 - sp - sm + g - te; dst[b + 8] = f8 - sp + sm - g - te; }
{ const sp = wp * (0.5 * (f9 + f10) - e4), sm = wm * (0.5 * (f9 - f10) - w2 * 3 * c4), g = kf * (Fx - Fy) * 1 / 36;
dst[b + 9] = f9 - sp - sm + g - te; dst[b + 10] = f10 - sp + sm - g - te; }
{ const sp = wp * (0.5 * (f11 + f12) - e5), sm = wm * (0.5 * (f11 - f12) - w2 * 3 * c5), g = kf * (Fx + Fz) * 1 / 36;
dst[b + 11] = f11 - sp - sm + g - te; dst[b + 12] = f12 - sp + sm - g - te; }
{ const sp = wp * (0.5 * (f13 + f14) - e6), sm = wm * (0.5 * (f13 - f14) - w2 * 3 * c6), g = kf * (Fx - Fz) * 1 / 36;
dst[b + 13] = f13 - sp - sm + g - te; dst[b + 14] = f14 - sp + sm - g - te; }
{ const sp = wp * (0.5 * (f15 + f16) - e7), sm = wm * (0.5 * (f15 - f16) - w2 * 3 * c7), g = kf * (Fy + Fz) * 1 / 36;
dst[b + 15] = f15 - sp - sm + g - te; dst[b + 16] = f16 - sp + sm - g - te; }
{ const sp = wp * (0.5 * (f17 + f18) - e8), sm = wm * (0.5 * (f17 - f18) - w2 * 3 * c8), g = kf * (Fy - Fz) * 1 / 36;
dst[b + 17] = f17 - sp - sm + g - te; dst[b + 18] = f18 - sp + sm - g - te; }
}
}
Mass through boundaries and force on walls
64-bit
Sums the net mass rate per boundary group and the momentum-exchange force on the part.
mass = Σ (f_in − f_out); F = −Σ cᵢ (f_in + f_out)
Running 64-bit sums over all boundary links in link order, without compensation. The mass sum adds terms of both signs, so the net through a group with little net flow has fewer correct digits than the gross flow.
flow.js
lines 404 to 417
PartCraft 0.8.1
/**
* Boundary bookkeeping from the current post-collision populations: per group the net mass rate into the fluid
* (lattice units, + = in), the number of links, and the force on wall links counted for forces (momentum exchange).
*/
export function boundaryStats(S) {
const { lat, src, ux, uy, uz } = S, L = lat.links, ng = lat.groups.length;
const mass = new Float64Array(ng), links = new Int32Array(ng), force = [0, 0, 0];
for (let k = 0; k < L.n; k++) {
const n = L.cell[k], q = L.dir[k], o = OPP[q], fin = incoming(lat, k, src, ux, uy, uz), fout = src[n * Q + o];
mass[L.g[k]] += L.kind[k] === K_BOUZIDI && lat.massFix !== false ? 0 : fin - fout; links[L.g[k]]++;
if (L.force[k]) { const s = fin + fout; force[0] -= CX[q] * s; force[1] -= CY[q] * s; force[2] -= CZ[q] * s; }
}
return { mass, links, force };
}
Stepping to a steady state
64-bit
Repeats time steps until the velocity field stops changing.
res = √(Σ |u − u_prev|² ÷ Σ |u|²) ÷ steps between checks
The residual is a relative change per step, checked every 100 steps against the tolerance (1e-7 by default). Reaching the step limit returns converged false; a sum that is not finite returns a blew-up flag instead of a result.
flow.js
lines 419 to 453
PartCraft 0.8.1
/**
* Time-step to a steady state. opts: { maxSteps, tol (relative velocity change per step), every (check interval),
* minSteps, onCheck({step, res}) → true to stop early (cancel), accept() → false to keep going although res < tol,
* ramp (steps over which inflow speeds rise from 0) }.
* Returns { steps, converged, res, history: [{step, res}], blewUp }.
*/
export function run(S, { maxSteps = 20000, tol = 1e-7, every = 100, minSteps = 0, onCheck = null, accept = null, ramp = 0 } = {}) {
const nf = S.lat.nf, prev = new Float64Array(nf * 3), hist = [];
let res = Infinity, converged = false;
const snap = () => { for (let n = 0; n < nf; n++) { prev[3 * n] = S.ux[n]; prev[3 * n + 1] = S.uy[n]; prev[3 * n + 2] = S.uz[n]; } };
// ramp: the inflow links' speeds rise smoothly (½ − ½cos) over `ramp` steps — a gentle start for runs near τ = ½
let rl = null, ra = null;
if (ramp > 0) { const L = S.lat.links, ids = []; for (let m = 0; m < S.lin.k.length; m++) if (L.kind[S.lin.k[m]] === K_VEL) ids.push(m); rl = Int32Array.from(ids); ra = Float64Array.from(ids, (m) => S.lin.a[m]); }
const stepR = () => {
if (rl && S.step <= ramp) { const f = 0.5 - 0.5 * Math.cos(Math.PI * Math.min(1, S.step / ramp)); for (let i = 0; i < rl.length; i++) S.lin.a[rl[i]] = f * ra[i]; }
step(S);
};
stepR(); snap();
while (S.step < maxSteps) {
const k = Math.min(every, maxSteps - S.step);
for (let s = 0; s < k; s++) stepR();
let d = 0, m = 0;
for (let n = 0; n < nf; n++) {
const a = S.ux[n], b = S.uy[n], c = S.uz[n];
d += (a - prev[3 * n]) ** 2 + (b - prev[3 * n + 1]) ** 2 + (c - prev[3 * n + 2]) ** 2; m += a * a + b * b + c * c;
}
if (!Number.isFinite(d) || !Number.isFinite(m)) return { steps: S.step, converged: false, res: NaN, history: hist, blewUp: true };
res = m > 0 ? Math.sqrt(d / m) / k : (d > 0 ? Infinity : 0);
hist.push({ step: S.step, res });
snap();
if (onCheck && onCheck({ step: S.step, res })) break;
if (S.step >= minSteps && res < tol && (!accept || accept())) { converged = true; break; }
}
return { steps: S.step, converged, res, history: hist, blewUp: false };
}
Lattice units from physical ones
64-bit
Chooses the lattice speed, viscosity, relaxation time and time step for a cell size, speed and fluid.
Re_cell = U dx ÷ ν; u = min(0.05, √(0.03 ÷ k), 0.5 Re_cell); τ = 3 u ÷ Re_cell + ½; dt = u dx ÷ U
64-bit closed form. The density budget uses the pipe friction factors 64 ÷ Re and, for turbulent runs, the Blasius fit 0.316 Re^−0.25; these set how fast the run is allowed to be, not the result's digits. A laminar case with τ under 0.5005 is refused.
flow.js
lines 457 to 476
PartCraft 0.8.1
/**
* Lattice parameters from physical ones. dx (m), U (m/s, the reference speed: inlet or free stream), nu (m²/s),
* Lc / Dc (flow path length and the narrowest width, in cells: they bound the density change), opts {uMax 0.05, nuMax 0.5, drho 0.03}.
* → { u, nu, tau, dt, Reg (cell Reynolds number), cu (m/s per lattice speed), cp (Pa per lattice density·ρ), estDrho, ok, why }
*/
export function scaling({ dx, U, nu, Lc = 50, Dc = 10, uMax = 0.05, nuMax = 0.5, drho = 0.03, internal = true, turbulent = false }) {
const Reg = U * dx / nu, Re = Reg * Dc;
// density change budget: pipe-friction loss over the path (Hagen–Poiseuille 64/Re, or Blasius 0.316 Re^−¼ when turbulent)
// + a few dynamic heads (internal); dynamic heads (external)
const f = turbulent ? Math.max(64 / Math.max(Re, 1e-12), 0.316 * Math.max(Re, 1) ** -0.25) : 64 / Math.max(Re, 1e-12);
const k = internal ? 3 * (f / 2 * Lc / Dc + 2) : 3 * 1.5;
let u = Math.min(uMax, Math.sqrt(drho / k), turbulent ? Infinity : nuMax * Reg);
const nuL = u / Reg, tau = 3 * nuL + 0.5, dt = u * dx / U;
const out = { u, nu: nuL, tau, dt, dx, Reg, Re, cu: dx / dt, estDrho: k * u * u, ok: true, why: '' };
// laminar: the cell must resolve the viscosity itself; turbulent: the eddy viscosity carries the stress, the molecular τ
// may sit at ½ (refused only when even a wall cell would be far outside the wall law's checked range)
if (!turbulent && tau < 0.5 + 5e-4) { out.ok = false; out.why = t('study.flow.err.cells_too_large', { re: Reg.toPrecision(3), n: Math.ceil(Reg / 30) }); }
if (turbulent && Reg > 1e5) { out.ok = false; out.why = t('study.flow.err.cells_too_large_turb', { re: Reg.toPrecision(3), n: Math.ceil(Reg / 2e4) }); }
return out;
}
In detail
Method
D3Q19 lattice Boltzmann with two-relaxation-time collision and Λ = 3/16; halfway or interpolated bounce-back walls (the file header cites Bouzidi, Firdaouss and Lallemand 2001, Ladd 1994, Ginzburg 2008, Izquierdo and Fueyo 2008); forces by momentum exchange.
Number types
Float64Array populations, density and velocity; Int8Array lattice directions; Int32Array indexes. Wall-link normals for the turbulence model are stored in a Float32Array.
Tolerances
Steady state: relative velocity change per step under 1e-7 by default (a study's default is 1e-6), checked every 100 steps. Interpolated wall used when the wall distance differs from ½ by more than 1e-6; the distance is clamped to [0.001, 1]. Lattice speed at most 0.05, estimated density change at most 3 percent, laminar relaxation time at least 0.5005.
Limits
Laminar, weakly compressible, steady results only (file header). Air: 1.2041 kg/m³ and 1.813e-5 Pa·s; water: 998.21 kg/m³ and 1.0016e-3 Pa·s.
Approach: 64-bit, stated limit
The arithmetic is 64-bit, but the method or the data sets the accuracy: a solver tolerance, a mesh, a step size, a datum shift, a table. The limit is stated where the result is shown.