Resources Numerical accuracy Calculation policy
Computes laminar flow and tracer mixing in tanks and flow-through boxes on a lattice, with a mixing time and a residence time distribution.
Distributions and fields are Float64Array and the arithmetic is 64-bit. The accuracy is set by the lattice method: a fixed cell and step of 1, a second-order equilibrium that the file header says needs speeds below 0.1, and sampling of results every few steps. The VTK export writes 6 significant digits.
Open item. Write the VTK values at full precision and declare them double.
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.
Tests tools/process-test/test-dec.mjs
The nine directions of the lattice and the weight of each.
w = 4/9, 1/9 (×4), 1/36 (×4)
The weights are 64-bit quotients of whole numbers, written as divisions and not as rounded decimals.
export const CX = [0, 1, 0, -1, 0, 1, -1, -1, 1], CY = [0, 0, 1, 0, -1, 1, 1, -1, -1], W = [4 / 9, 1 / 9, 1 / 9, 1 / 9, 1 / 9, 1 / 36, 1 / 36, 1 / 36, 1 / 36], OPP = [0, 3, 4, 1, 2, 7, 8, 5, 6];Creates the grid, the solid mask, the force field and the starting distributions.
The grid counts are made whole with | 0; they are counts, not measurements. Every field is a Float64Array. Both relaxation times must be above ½, and the grid is limited to four million cells.
export function createMixer(o) {
const nx = o.nx | 0, ny = o.ny | 0; if (!(nx >= 4 && ny >= 4 && nx * ny <= 4e6)) throw new Error('the mixer grid is 4 x 4 up to four million cells');
const N = nx * ny, tau = o.tau ?? 0.8, tauD = o.tauD ?? 0.7; if (!(tau > 0.5) || !(tauD > 0.5)) throw new Error('tau and tauD must exceed 1/2');
const solid = new Uint8Array(N);
if (typeof o.solid === 'function') { for (let y = 0; y < ny; y++) for (let x = 0; x < nx; x++) solid[x + nx * y] = o.solid(x, y) ? 1 : 0; } else if (o.solid) solid.set(o.solid);
const S = { nx, ny, N, tau, tauD, solid, f: new Float64Array(9 * N), g: new Float64Array(9 * N), f2: new Float64Array(9 * N), g2: new Float64Array(9 * N), rho: new Float64Array(N).fill(1), ux: new Float64Array(N), uy: new Float64Array(N), C: new Float64Array(N),
force: o.force || [0, 0], sources: (o.sources || []).map((s) => ({ ...s })), flow: o.flow || null, t: 0 };
for (let i = 0; i < N; i++) for (let k = 0; k < 9; k++) { S.f[9 * i + k] = W[k]; S.g[9 * i + k] = 0; }
// force per cell (Guo): uniform + the sources' discs
S.Fx = new Float64Array(N).fill(S.force[0]); S.Fy = new Float64Array(N).fill(S.force[1]);
for (const s of S.sources) for (let y = 0; y < ny; y++) for (let x = 0; x < nx; x++) if (!solid[x + nx * y] && (x - s.x) ** 2 + (y - s.y) ** 2 <= s.r * s.r) { S.Fx[x + nx * y] += s.fx || 0; S.Fy[x + nx * y] += s.fy || 0; }
if (S.flow) { const u = S.flow.u; for (let i = 0; i < N; i++) { if (solid[i]) continue; for (let k = 0; k < 9; k++) S.f[9 * i + k] = W[k] * (1 + 3 * CX[k] * u + 3 * u * u * (CX[k] * CX[k] - 1 / 3) * 1.5); S.ux[i] = u; } }
return S;
}The equilibrium the flow and the tracer relax toward.
f_eq = w ρ (1 + 3 c·u + 4.5 (c·u)² − 1.5 u²); g_eq = w C (1 + 3 c·u)
One 64-bit expression each. The flow equilibrium is second order in velocity and the tracer's is first order, which is the method's own truncation.
const feq = (k, rho, ux, uy) => { const cu = CX[k] * ux + CY[k] * uy; return W[k] * rho * (1 + 3 * cu + 4.5 * cu * cu - 1.5 * (ux * ux + uy * uy)); };
const geq = (k, C, ux, uy) => W[k] * C * (1 + 3 * (CX[k] * ux + CY[k] * uy));Density, velocity and tracer value of each cell, and the relaxation of its distributions.
ρ = Σ f; u = (Σ f c + ½ F) ÷ ρ; f ← f − (f − f_eq) ÷ τ + source
Density and momentum are sums of the nine populations in a fixed order, and the velocity adds half the force (Guo). The division by density is not guarded; density starts at 1 in every fluid cell.
let r = 0, mx = 0, my = 0, c = 0; const b = 9 * i;
for (let k = 0; k < 9; k++) { const fk = S.f[b + k]; r += fk; mx += fk * CX[k]; my += fk * CY[k]; c += S.g[b + k]; }
const ux = (mx + 0.5 * Fx[i]) / r, uy = (my + 0.5 * Fy[i]) / r; S.rho[i] = r; S.ux[i] = ux; S.uy[i] = uy; S.C[i] = c;
for (let k = 0; k < 9; k++) {
const cu = CX[k] * ux + CY[k] * uy, Sk = (1 - 0.5 * omega) * W[k] * (3 * ((CX[k] - ux) * Fx[i] + (CY[k] - uy) * Fy[i]) + 9 * cu * (CX[k] * Fx[i] + CY[k] * Fy[i]));
S.f[b + k] += -omega * (S.f[b + k] - feq(k, r, ux, uy)) + Sk;
S.g[b + k] += -omegaD * (S.g[b + k] - geq(k, c, ux, uy));Sets the last column of a flow-through box from its neighbor.
f_out = f_eq(1, u_nb) + (f_nb − f_eq(ρ_nb, u_nb))
The outlet takes the equilibrium at density 1 and the neighbor's velocity, plus the neighbor's non-equilibrium part, which is a difference of two close numbers. The tracer is copied from the neighbor.
let r = 0, mx = 0, my = 0; for (let k = 0; k < 9; k++) { const v = S.f2[9 * iM + k]; r += v; mx += v * CX[k]; my += v * CY[k]; }
const ux = mx / r, uy = my / r;
for (let k = 0; k < 9; k++) { S.f2[9 * iN + k] = feq(k, 1, ux, uy) + (S.f2[9 * iM + k] - feq(k, r, ux, uy)); S.g2[9 * iN + k] = S.g2[9 * iM + k]; }The velocity reported at a position in the grid.
Returns the values of the nearest cell (Math.round of x and y). Nothing is interpolated, so the result changes in steps of one cell.
export const velocityAt = (S, x, y) => { const i = Math.round(x) + S.nx * Math.round(y); return [S.ux[i], S.uy[i]]; };A copy of the velocity, speed, density or tracer field for reading or plotting.
speed = √(u_x² + u_y²)
The copy is a Float64Array with solid cells set to 0. Speed is Math.hypot of the two components.
export function field(S, name) {
if (name === 'speed') return Float64Array.from(S.ux, (v, i) => (S.solid[i] ? 0 : Math.hypot(v, S.uy[i])));
const a = S[name]; if (!a) throw new Error('field: ux | uy | C | rho | speed'); return Float64Array.from(a, (v, i) => (S.solid[i] ? 0 : v));
}The number of steps until the tracer is evenly spread.
deviation = max |C − mean|; reached when deviation ≤ level × starting deviation
The mean is a sum over the fluid cells in order and the deviation is the largest absolute difference from it. It is tested every 50 steps by default, so the step count returned is that of the first check at or under the level.
const stats = (S) => { let sum = 0, n = 0; for (let i = 0; i < S.N; i++) if (!S.solid[i]) { sum += S.C[i]; n++; } const mean = sum / n; let dev = 0; for (let i = 0; i < S.N; i++) if (!S.solid[i]) dev = Math.max(dev, Math.abs(S.C[i] - mean)); return { mean, dev }; };
export function mixingTime(S, { level = 0.05, maxSteps = 20000, every = 50 } = {}) {
step(S, 0); const s0 = stats(S), history = [[S.t, s0.dev]]; if (!(s0.dev > 0)) return { steps: 0, reached: true, history };
const t0 = S.t;
for (let n = 0; n < maxSteps; n += every) {
step(S, every); const s = stats(S); history.push([S.t - t0, s.dev]);
if (s.dev <= level * s0.dev) return { steps: S.t - t0, reached: true, history, mean: s.mean };
}
return { steps: null, reached: false, history };
}The distribution of times a pulse of tracer takes to leave a flow-through box, with its mean, variance and the nominal time.
E = flux ÷ Σ(flux × Δt); mean = Σ t E Δt; variance = Σ t² E Δt − mean²; nominal = V ÷ Q
The outlet flux Σ C u_x in the next-to-last column is sampled every few steps and summed as rectangles. The variance is a difference of two larger sums. The nominal time is the count of fluid cells divided by the outlet flow.
/** Residence time distribution of a flow-through box: a one-step pulse of tracer at the inlet, then the tracer flux through the last column. E(t) is normalised to unit area. */
export function rtd(S, { steps = 4000, every = 5, pulse = 1 } = {}) {
if (!S.flow) throw new Error('rtd needs a flow-through mixer (flow: { u })');
const t = [], E = [], c0 = S.flow.C; let area = 0;
S.flow.C = pulse; step(S, 1); S.flow.C = 0;
const base = S.t;
for (let n = 0; n < steps; n += every) {
step(S, every); let flux = 0; for (let y = 0; y < S.ny; y++) { const i = S.nx - 2 + S.nx * y; if (!S.solid[i]) flux += S.C[i] * S.ux[i]; }
t.push(S.t - base); E.push(flux); area += flux * every;
}
S.flow.C = c0;
const En = E.map((v) => v / (area || 1)); let mean = 0, m2 = 0; for (let k = 0; k < t.length; k++) { mean += t[k] * En[k] * every; m2 += t[k] * t[k] * En[k] * every; }
let Q = 0, V = 0; for (let y = 0; y < S.ny; y++) { if (!S.solid[S.nx - 2 + S.nx * y]) Q += S.ux[S.nx - 2 + S.nx * y]; for (let x = 1; x <= S.nx - 2; x++) if (!S.solid[x + S.nx * y]) V++; }
return { t, E: En, mean, variance: m2 - mean * mean, area, nominal: Q > 0 ? V / Q : null }; // nominal = V/Q from the outlet flux: the mean of a complete RTD equals it
}Writes the velocity and tracer fields as a legacy ASCII VTK file.
Every value is written with toPrecision(6), so the file carries 6 significant digits of the 64-bit result. The header declares the arrays as float.
export function vtk(S, name = 'mixer') {
const L = ['# vtk DataFile Version 3.0', name, 'ASCII', 'DATASET STRUCTURED_POINTS', `DIMENSIONS ${S.nx} ${S.ny} 1`, 'ORIGIN 0 0 0', 'SPACING 1 1 1', `POINT_DATA ${S.N}`, 'VECTORS velocity float'];
for (let i = 0; i < S.N; i++) L.push(`${(S.solid[i] ? 0 : S.ux[i]).toPrecision(6)} ${(S.solid[i] ? 0 : S.uy[i]).toPrecision(6)} 0`);
L.push('SCALARS tracer float 1', 'LOOKUP_TABLE default'); for (let i = 0; i < S.N; i++) L.push((S.solid[i] ? 0 : S.C[i]).toPrecision(6)); return L.join('\n') + '\n';
}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.