Computes laminar flow in a rectangular region and the mixing of a tracer in it: velocities, pressure, concentration and a mixing time.
How the precision is kept
All fields are Float64Array and every operation is 64-bit. The accuracy is set by the method: a uniform grid, an explicit time step chosen from stability limits, a hybrid central and upwind scheme, and a pressure solve that stops at a relative residual of 1e-8. The file header states it is for laminar flow only.
Open item. Declare the VTK arrays as double so a reader keeps the 64-bit values that are written.
Calculations under this policy
14 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.
The value used just outside the region so that the velocity at a wall equals the wall's own.
u_ghost = 2 × u_wall − u_inside
The ghost value mirrors the first inner value about the wall velocity, so their mean at the wall is the wall velocity. At an inflow side the tangential value is mirrored to zero, and at an outflow side it is copied.
cfd/fv.js
line 41
CadCraft 0.8.2
if (j < 0) { const b = bc.bottom; if (b.type === 'wall') return 2 * (b.u || 0) - U[i + (nx + 1) * 0]; if (b.type === 'inflow') return -U[i]; return U[i]; }
Value at a cell face
64-bit
The value carried across a cell face by the flow.
central ½(a + b) when |u| h ÷ ν ≤ 2, else the upwind value
The switch is a hard threshold at a cell Peclet number of 2, so the scheme drops to first-order upwind wherever the flow is fast for the grid.
cfd/fv.js
lines 51 to 52
CadCraft 0.8.2
/** hybrid face value: central unless the cell Peclet number is above 2, then upwind */
face(a, b, vel, h) { return Math.abs(vel) * h / this.nu > 2 ? (vel > 0 ? a : b) : 0.5 * (a + b); }
Time step
64-bit
The largest stable time step for the current velocities.
u_max is the largest absolute face velocity, floored at 1e-12, and h the smaller cell side. The three limits are the advective, the diffusive and the explicit central-scheme limits named in the code's comment.
cfd/fv.js
lines 53 to 56
CadCraft 0.8.2
dtLimit() {
let um = 1e-12; for (const a of this.U) um = Math.max(um, Math.abs(a)); for (const a of this.V) um = Math.max(um, Math.abs(a));
const h = Math.min(this.hx, this.hy); return Math.min(0.4 * h / um, 0.2 * h * h / this.nu, 0.8 * this.nu / (um * um)); // advective CFL, diffusion limit, and the explicit-central stability limit dt < 2 nu / u²
}
Momentum step for u
64-bit
The predicted x velocity at each face after one time step.
u* = u + dt (−convection + ν ∇²u + source)
Explicit Euler in time with second-order central differences for diffusion and the hybrid face values for convection. Values are read from Float64Array and combined in 64-bit; the error is that of the grid and the time step.
Sets up the pressure equation that makes the velocity field divergence-free.
∇²p = (∇·u*) ÷ dt
The Laplacian is applied without a stored matrix, with p = 0 half a cell beyond outflow sides and zero gradient elsewhere. Without an outflow side the mean of the right side is subtracted so the equation has a solution.
cfd/fv.js
lines 87 to 100
CadCraft 0.8.2
for (let j = 0; j < ny; j++) for (let i = 0; i < nx; i++) b[i + nx * j] = ((U[i + 1 + (nx + 1) * j] - U[i + (nx + 1) * j]) / hx + (V[i + nx * (j + 1)] - V[i + nx * j]) / hy) / dt;
// matrix-free Laplacian; Dirichlet p = 0 on outflow sides (distance h/2), Neumann elsewhere
const bc = this.bc, dirL = bc.left.type === 'outflow', dirR = bc.right.type === 'outflow', dirB = bc.bottom.type === 'outflow', dirT = bc.top.type === 'outflow';
const A = (x, y) => {
for (let j = 0; j < ny; j++) for (let i = 0; i < nx; i++) {
const k = i + nx * j, c = x[k]; let s = 0;
if (i > 0) s += (x[k - 1] - c) / (hx * hx); else if (dirL) s += (-2 * c) / (hx * hx);
if (i < nx - 1) s += (x[k + 1] - c) / (hx * hx); else if (dirR) s += (-2 * c) / (hx * hx);
if (j > 0) s += (x[k - nx] - c) / (hy * hy); else if (dirB) s += (-2 * c) / (hy * hy);
if (j < ny - 1) s += (x[k + nx] - c) / (hy * hy); else if (dirT) s += (-2 * c) / (hy * hy);
y[k] = -s; // positive semi-definite form
}
};
const rhs = b.map((v) => -v); if (!this.outflow) { let m = 0; for (const v of rhs) m += v; m /= n; for (let k = 0; k < n; k++) rhs[k] -= m; }
Pressure solve
64-bit
Solves the pressure equation by conjugate gradient.
Unpreconditioned conjugate gradient starting from the previous pressure. It stops when the squared residual falls below 1e-16 of its starting value (a relative residual of 1e-8), after 4(nx + ny) + 200 passes, or when pᵀAp is not positive. No flag is returned when the pass limit ends it.
cfd/fv.js
lines 102 to 108
CadCraft 0.8.2
const x = P, r = new Float64Array(n), p = new Float64Array(n), Ap = new Float64Array(n); A(x, Ap); for (let k = 0; k < n; k++) { r[k] = rhs[k] - Ap[k]; p[k] = r[k]; }
let rr = 0; for (const v of r) rr += v * v; const rr0 = rr || 1;
for (let it = 0; it < 4 * (nx + ny) + 200 && rr > 1e-16 * rr0; it++) {
A(p, Ap); let pAp = 0; for (let k = 0; k < n; k++) pAp += p[k] * Ap[k]; if (pAp <= 0) break;
const al = rr / pAp; let rn = 0; for (let k = 0; k < n; k++) { x[k] += al * p[k]; r[k] -= al * Ap[k]; rn += r[k] * r[k]; }
const be = rn / rr; rr = rn; for (let k = 0; k < n; k++) p[k] = r[k] + be * p[k];
}
Velocity correction
64-bit
Removes the pressure gradient from the predicted velocities.
u = u* − dt (p_right − p_left) ÷ h
One subtraction per face. At an outflow side the gradient is taken over half a cell to p = 0.
Fluxes use the van Leer limiter, with r set to 0 when the downwind difference is exactly 0, plus central diffusion. The update is in flux form: what leaves one cell through an inner face enters its neighbor.
cfd/fv.js
lines 123 to 136
CadCraft 0.8.2
advectScalar(s, dt) {
const { nx, ny, hx, hy, U, V } = this, C = s.C, N = C.slice(), bc = this.bc;
const lim = (r) => (r + Math.abs(r)) / (1 + Math.abs(r)); // van Leer
const cAt = (i, j) => (i < 0 ? (bc.left.type === 'inflow' ? bc.left.c?.(this.t) ?? 0 : C[nx * j]) : i >= nx ? C[nx - 1 + nx * j] : j < 0 ? C[i] : j >= ny ? C[i + nx * (ny - 1)] : C[i + nx * j]);
const fx = (i, j) => { const u = U[i + (nx + 1) * j]; if (i === 0 && bc.left.type !== 'outflow') return u * (bc.left.type === 'inflow' ? (bc.left.c?.(this.t) ?? 0) : 0); const a = cAt(i - 1, j), b = cAt(i, j); const up = u > 0 ? a : b, dn = u > 0 ? b : a, far = u > 0 ? cAt(i - 2, j) : cAt(i + 1, j); const r = dn - up === 0 ? 0 : (up - far) / (dn - up); return u * (up + 0.5 * lim(r) * (dn - up)) - s.D * (b - a) / hx; };
const fy = (i, j) => { const v = V[i + nx * j]; if (j === 0 || j === ny) { if (j === 0 && bc.bottom.type === 'inflow') return v * (bc.bottom.c?.(this.t) ?? 0); if (j === ny && bc.top.type === 'inflow') return v * (bc.top.c?.(this.t) ?? 0); if ((j === 0 && bc.bottom.type !== 'outflow') || (j === ny && bc.top.type !== 'outflow')) return 0; } const a = cAt(i, j - 1), b = cAt(i, j); const up = v > 0 ? a : b, dn = v > 0 ? b : a, far = v > 0 ? cAt(i, j - 2) : cAt(i, j + 1); const r = dn - up === 0 ? 0 : (up - far) / (dn - up); return v * (up + 0.5 * lim(r) * (dn - up)) - s.D * (b - a) / hy; };
for (let j = 0; j < ny; j++) for (let i = 0; i < nx; i++) {
const wall = (side) => side.type === 'wall';
const flxW = (i === 0 && wall(bc.left)) ? 0 : fx(i, j), flxE = (i === nx - 1 && wall(bc.right)) ? 0 : (i === nx - 1 ? U[nx + (nx + 1) * j] * C[i + nx * j] : fx(i + 1, j));
const flyS = (j === 0 && wall(bc.bottom)) ? 0 : fy(i, j), flyN = (j === ny - 1 && wall(bc.top)) ? 0 : (j === ny - 1 ? V[i + nx * ny] * C[i + nx * j] : fy(i, j + 1));
N[i + nx * j] = C[i + nx * j] - dt * ((flxE - flxW) / hx + (flyN - flyS) / hy) + (s.source ? dt * s.source((i + 0.5) * hx, (j + 0.5) * hy, this.t) : 0);
}
C.set(N);
}
Run to a steady state
64-bit
Steps the flow until it stops changing or a limit is reached.
change = max |u′ − u| ÷ dt
The run stops when the largest change of any velocity component per unit time is below the tolerance (1e-7 by default, absolute) after more than 10 steps, or at the step or time limit. The last change is returned.
cfd/fv.js
lines 137 to 146
CadCraft 0.8.2
/** Run until the velocity change per unit time is under `steady`, or maxSteps / tEnd. → { steps, t, change } */
run({ steady = 1e-7, maxSteps = 100000, tEnd = Infinity } = {}) {
let ch = Infinity;
while (this.steps < maxSteps && this.t < tEnd) {
const dt = Math.min(this.dtLimit(), tEnd - this.t), u0 = this.U.slice(), v0 = this.V.slice(); this.step(dt);
ch = 0; for (let k = 0; k < u0.length; k++) ch = Math.max(ch, Math.abs(this.U[k] - u0[k])); for (let k = 0; k < v0.length; k++) ch = Math.max(ch, Math.abs(this.V[k] - v0[k])); ch /= dt;
if (ch < steady && this.steps > 10) break;
}
return { steps: this.steps, t: this.t, change: ch };
}
Each component is the mean of its two face values, stored in a Float64Array.
cfd/fv.js
lines 147 to 148
CadCraft 0.8.2
/** u, v at cell centres */
cellVelocity() { const { nx, ny, U, V } = this, u = new Float64Array(nx * ny), v = new Float64Array(nx * ny); for (let j = 0; j < ny; j++) for (let i = 0; i < nx; i++) { u[i + nx * j] = 0.5 * (U[i + (nx + 1) * j] + U[i + 1 + (nx + 1) * j]); v[i + nx * j] = 0.5 * (V[i + nx * j] + V[i + nx * (j + 1)]); } return { u, v }; }
Velocity at a point
64-bit
The velocity at any point of the region.
Bilinear interpolation of each staggered component. The fractional index is clamped to the grid, less 1e-12 at the upper end, so points outside take the nearest edge values.
cfd/fv.js
lines 149 to 153
CadCraft 0.8.2
/** Velocity at a point (bilinear on the staggered components) */
sample(x, y) {
const { nx, ny, hx, hy } = this, bil = (arr, w, h, ox, oy, X, Y) => { let fx = X / hx - ox, fy = Y / hy - oy; fx = Math.min(Math.max(fx, 0), w - 1 - 1e-12); fy = Math.min(Math.max(fy, 0), h - 1 - 1e-12); const i = Math.floor(fx), j = Math.floor(fy), a = fx - i, b = fy - j; return (1 - a) * (1 - b) * arr[i + w * j] + a * (1 - b) * arr[i + 1 + w * j] + (1 - a) * b * arr[i + w * (j + 1)] + a * b * arr[i + 1 + w * (j + 1)]; };
return [bil(this.U, nx + 1, ny, 0, 0.5, x, y), bil(this.V, nx, ny + 1, 0.5, 0, x, y)];
}
Stream function
64-bit
The stream function at the cell corners.
ψ(i, j) = ψ(i, j − 1) + U(i, j − 1) × h_y
A running sum up each column from ψ = 0 at the bottom, so rounding accumulates with the number of rows.
cfd/fv.js
lines 154 to 155
CadCraft 0.8.2
/** Stream function psi at the lower-left cell corners by integrating U upward: psi(i, j) with psi = 0 on the bottom */
streamFunction() { const { nx, ny, hy, U } = this, psi = new Float64Array((nx + 1) * (ny + 1)); for (let i = 0; i <= nx; i++) for (let j = 1; j <= ny; j++) psi[i + (nx + 1) * j] = psi[i + (nx + 1) * (j - 1)] + U[i + (nx + 1) * (j - 1)] * hy; return psi; }
VTK export
64-bit
Writes velocity, pressure and tracer fields as a legacy ASCII VTK file.
Values are written with String(), the shortest text that reads back to the same 64-bit number. The file header declares the arrays as float, the 32-bit type of the legacy VTK format, so a reader that follows the header stores them in 32-bit.
cfd/fv.js
lines 156 to 163
CadCraft 0.8.2
/** Legacy ASCII VTK (an open format) of cell data: velocity, pressure, scalars */
vtk(title = 'cadcraft flow') {
const { nx, ny, hx, hy } = this, { u, v } = this.cellVelocity(); const L = ['# vtk DataFile Version 3.0', title, 'ASCII', 'DATASET STRUCTURED_POINTS', `DIMENSIONS ${nx + 1} ${ny + 1} 1`, `ORIGIN 0 0 0`, `SPACING ${hx} ${hy} 1`, `CELL_DATA ${nx * ny}`, 'VECTORS velocity float'];
for (let k = 0; k < nx * ny; k++) L.push(`${u[k]} ${v[k]} 0`);
L.push('SCALARS pressure float 1', 'LOOKUP_TABLE default'); for (let k = 0; k < nx * ny; k++) L.push(String(this.P[k]));
for (const s of Object.values(this.scalars)) { L.push(`SCALARS ${s.name} float 1`, 'LOOKUP_TABLE default'); for (let k = 0; k < nx * ny; k++) L.push(String(s.C[k])); }
return L.join('\n') + '\n';
}
Mixing time
64-bit
The time at which a tracer's spread falls below a target.
CV = √(Σ (c − mean)² ÷ n) ÷ mean
The mean and the variance are sums over all cells in order. The spread is tested every 20 steps by default, so the time returned is that of the first check at or below the target, not the exact crossing.
cfd/fv.js
lines 166 to 172
CadCraft 0.8.2
/** Mixing statistics of a scalar: mean, coefficient of variation, and the first time the CV falls below `target` while stepping. */
export function mixingTime(flow, name, { target = 0.05, tMax = 1e9, check = 20 } = {}) {
const s = flow.scalars[name], cv = () => { let m = 0; for (const c of s.C) m += c; m /= s.C.length; let v = 0; for (const c of s.C) v += (c - m) ** 2; return m > 0 ? Math.sqrt(v / s.C.length) / m : Infinity; };
const hist = [[flow.t, cv()]];
while (flow.t < tMax && hist[hist.length - 1][1] > target) { for (let k = 0; k < check; k++) flow.step(); hist.push([flow.t, cv()]); }
return { time: hist[hist.length - 1][1] <= target ? flow.t : null, history: hist };
}
In detail
Method
Staggered grid with a fractional-step projection; hybrid differencing by cell Peclet number; conjugate gradient for the pressure Poisson equation; van Leer limited advection for the tracer (the file header cites Harlow and Welch 1965, Chorin 1968, Spalding 1972, Sweby 1984).
Number types
Float64Array for velocity, pressure and tracer fields; 64-bit scalars.
Tolerances
Pressure solve: squared residual below 1e-16 of its starting value, at most 4(nx + ny) + 200 passes. Steady state: largest velocity change per unit time below 1e-7 by default. Time step: min(0.4 h ÷ u, 0.2 h² ÷ ν, 0.8 ν ÷ u²). Hybrid switch at a cell Peclet number of 2.
Limits
Explicit first-order time stepping; first-order upwind where the Peclet number exceeds 2. The pressure solve does not report when it ends at the pass limit. The VTK export declares its arrays as float.
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.