Solves a looped water network for the head and pressure at every node and the flow, velocity and head loss in every pipe, then runs fire-flow and extended-period cases on it.
How the precision is kept
Everything is 64-bit: plain numbers and Float64Array. The solver works in SI and converts project units in and out with fixed factors; the gallon-per-minute and psi-per-foot factors are written to six significant digits. Head loss is the empirical Hazen-Williams formula, and the network is solved by an iteration that stops on a relative flow change of 1e-7, so the method and its constants set the accuracy, not the number type.
Open item. Write the gallon-per-minute factor as 3.785411784e-3 ÷ 60 and derive psi per foot from exact constants, and export the EPANET-format numbers at full precision instead of 4 decimals.
Calculations under this policy
12 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 factors that turn project units (ft, in, gpm, psi or m, mm, L/s, kPa) into the SI values the solver works in.
ft = 0.3048 m; in = 0.0254 m; gpm = 6.30902e-5 m³/s; psi per ft = 0.433527; kPa per m = 9.80665
0.3048, 0.0254 and g = 9.80665 are the exact defined values. GPM is written as 6.30902e-5, six significant digits; the exact value is 3.785411784e-3 ÷ 60 = 6.30901964e-5. PSI_PER_FT = 0.433527 is also a six-digit constant, and both pressure factors assume water at 1000 kg/m³.
water.js
lines 28 to 33
CivilCraft 0.8.1
const G = 9.80665;
const FT = 0.3048, INCH = 0.0254, GPM = 6.30902e-5, LS = 1e-3, MM = 1e-3;
export const PSI_PER_FT = 0.433527, KPA_PER_M = 9.80665;
export const LIMITS = { nodes: 5000, pipes: 8000, text: 120, pattern: 1000, curve: 3 };
const si = (units) => (units === 'm' ? { len: 1, dia: MM, q: LS, p: KPA_PER_M } : { len: FT, dia: INCH, q: GPM, p: PSI_PER_FT / FT });
Pump curve fit
64-bit
Turns one or three pump curve points into the head a pump adds at any flow.
H = H₀ − r·Qⁿ; n = ln((H₀ − H₂) ÷ (H₀ − H₁)) ÷ ln(Q₂ ÷ Q₁); r = (H₀ − H₁) ÷ Q₁ⁿ
Closed form in 64-bit with two logarithms and one power. A single point uses a shutoff head of 1.33334 × H₁ and zero head at 2 × Q₁. With three points whose first flow is not zero, the first point's head is used as the shutoff head at zero flow, so the fitted curve does not pass through that first point.
water.js
lines 91 to 116
CivilCraft 0.8.1
/** A pump curve's head at flow Q (SI): one point → shutoff 1.33 H and 2 Q at 0 head (the usual convention); three points →
* the power law H = H0 − r·Q^n through the shutoff point and the other two. Returns { h0, r, n }. */
export function pumpCurve(curve) {
const c = [...(curve || [])].sort((x, y) => x[0] - y[0]);
if (c.length === 1) {
const [q1, h1] = c[0];
if (!(q1 > 0 && h1 > 0)) return null;
const h0 = 1.33334 * h1, q2 = 2 * q1, h2 = 0;
const n = Math.log((h0 - h2) / (h0 - h1)) / Math.log(q2 / q1);
return { h0, r: (h0 - h1) / Math.pow(q1, n), n };
}
if (c.length === 3 && c[0][0] === 0) {
const [, h0] = c[0], [q1, h1] = c[1], [q2, h2] = c[2];
if (!(h0 > h1 && h1 > h2 && q2 > q1 && q1 > 0)) return null;
const n = Math.log((h0 - h2) / (h0 - h1)) / Math.log(q2 / q1);
return { h0, r: (h0 - h1) / Math.pow(q1, n), n };
}
if (c.length === 3) {
const [q0, h0] = c[0], [q1, h1] = c[1], [q2, h2] = c[2]; // a design point in the middle, as three points of a curve
if (!(h0 > h1 && h1 > h2 && q2 > q1 && q1 > q0)) return null;
const hh = h0 + (h0 - h1) * (q0 / Math.max(1e-12, q1 - q0)) * 0; // no shutoff given: use the first point as shutoff
const n = Math.log((hh - h2) / (hh - h1)) / Math.log(q2 / q1);
return { h0: hh, r: (hh - h1) / Math.pow(q1, n), n };
}
return null;
}
Hazen-Williams head loss
64-bit
The friction head loss in a pipe for a given flow, length, diameter and C factor, in SI.
h = 10.67 × L × |Q|^1.852 ÷ (C^1.852 × D^4.871) × sign(Q)
One 64-bit expression with three Math.pow calls. The coefficient 10.67 and the exponents 1.852 and 4.871 are the formula's own rounded empirical constants, so they limit the result before the arithmetic does.
water.js
lines 118 to 119
CivilCraft 0.8.1
/** Hazen–Williams in SI: head loss (m) for flow Q (m³/s) in a pipe L m long, D m wide. */
export const hwLoss = (Q, L, D, C) => 10.67 * L * Math.pow(Math.abs(Q), 1.852) / (Math.pow(C, 1.852) * Math.pow(D, 4.871)) * Math.sign(Q);
Pipe resistance and minor-loss coefficients
64-bit
The two coefficients each pipe carries into the solver: friction resistance and minor loss.
r = 10.67 × L ÷ (C^1.852 × D^4.871); m = 8K ÷ (π² g D⁴)
Computed once per pipe in SI from the length and diameter after unit conversion. A pipe with no stored length uses the plan distance between its nodes from Math.hypot, so elevation difference is not included in that length. π is Math.PI and g is 9.80665.
The head loss across a pipe, valve or pump at the current flow, and its derivative, for each solver pass.
pipe: f = (r·|q|^0.852 + m·|q|)·q, g = 1.852·r·|q|^0.852 + 2·m·|q|; pump: f = −(H₀ − r·qⁿ), g = r·n·qⁿ⁻¹
The flow used inside the powers is floored at 1e-7 m³/s so the gradient never reaches zero; the gradient is also floored at 1e-9 for pipes and 1e-6 for pumps. A closed link returns a loss of q × 1e8 with gradient 1e8.
water.js
lines 167 to 176
CivilCraft 0.8.1
const link = (k, q) => { // → [f(q) head loss a→b, g = df/dq]
if (k.status === 'closed') return [q * 1e8, 1e8];
if (k.kind === 'pump') {
const c = k.pc, qq = Math.max(q, QMIN);
return [-(c.h0 - c.r * Math.pow(qq, c.n)), Math.max(c.r * c.n * Math.pow(qq, c.n - 1), 1e-6)];
}
const aq = Math.max(Math.abs(q), QMIN);
const f = (k.r * Math.pow(aq, 0.852) + k.m * aq) * q;
return [f, Math.max(1.852 * k.r * Math.pow(aq, 0.852) + 2 * k.m * aq, 1e-9)];
};
Node equations for one pass
64-bit
Builds the linear system for the unknown heads from the link coefficients, demands and fixed heads, and solves it.
Q = P·(Hₐ − H_b) + Y; P = 1 ÷ g; Y = Q − f ÷ g; Σ flows at each node = demand
The diagonal and right-hand side are Float64Array and are summed link by link in input order. A closed link is still added with P = 1e-8, so it contributes a very small conductance to the matrix although its flow is then set to zero.
water.js
lines 201 to 214
CivilCraft 0.8.1
const diag = new Float64Array(M), rhs = new Float64Array(M), rows = Array.from({ length: M }, () => []);
const extra = new Float64Array(N);
for (const k of L) if (k.kind === 'prv' && k.status === 'active') { const qp = qprv.get(k.id) || 0; extra[k.a] += qp; }
for (let j = 0; j < M; j++) rhs[j] = -(dem[unk[j]] + extra[unk[j]]);
L.forEach((k, i) => {
if (!(P[i] > 0)) return;
const ua = pos0[k.a], ub = pos0[k.b];
// node b gains +Q, node a loses Q: Q = P (Ha − Hb) + Y
if (ub >= 0) { diag[ub] += P[i]; rhs[ub] += Y[i]; if (ua >= 0) rows[ub].push([ua, -P[i]]); else if (fh[k.a] != null) rhs[ub] += P[i] * fh[k.a]; }
if (ua >= 0) { diag[ua] += P[i]; rhs[ua] -= Y[i]; if (ub >= 0) rows[ua].push([ub, -P[i]]); else if (fh[k.b] != null) rhs[ua] += P[i] * fh[k.b]; }
});
const Hn = solveSpd(M, diag, rows, rhs, unk.map((i) => H[i]));
for (let j = 0; j < M; j++) H[unk[j]] = Hn[j];
for (let i = 0; i < N; i++) if (fh[i] != null) H[i] = fh[i];
Flow update and convergence test
64-bit
Updates every link flow from the new heads and decides when the network is solved.
stop when Σ|Q_new − Q_old| ÷ Σ|Q_new| < 1e−7 and no valve or pump changed status
The test is relative to the total flow, so it scales with the network. The loop runs at most 300 passes; if it has not met the test the result carries a noconverge problem and ok is false. Reported flows are therefore good to the 1e-7 relative change of the last pass, not to the full 64-bit width.
water.js
lines 216 to 231
CivilCraft 0.8.1
let dq = 0, sq = 0;
L.forEach((k, i) => {
let qn;
if (k.kind === 'prv' && k.status === 'active') qn = null;
else if (!reach[k.a] && !reach[k.b]) qn = 0;
else if (k.status === 'closed') qn = 0;
else qn = P[i] * (H[k.a] - H[k.b]) + Y[i];
if (qn != null) { dq += Math.abs(qn - Q[i]); sq += Math.abs(qn); Q[i] = qn; }
});
// active PRV flow = what continuity at its downstream node needs
for (const k of L) if (k.kind === 'prv' && k.status === 'active') {
let qp = dem[k.b];
L.forEach((o, i) => { if (o === k) return; if (o.a === k.b) qp += Q[i]; if (o.b === k.b) qp -= Q[i]; });
dq += Math.abs(qp - (Q[L.indexOf(k)] || 0)); sq += Math.abs(qp);
Q[L.indexOf(k)] = qp; qprv.set(k.id, qp);
}
Linear solve for the heads
64-bit
Solves the symmetric positive definite system for the unknown node heads.
A = L·Lᵀ, forward and back substitution; or conjugate gradient with z = r ÷ diag(A)
Up to 250 unknowns a dense in-place Cholesky on Float64Array rows is used; a pivot that comes out at or below zero is replaced by √1e-30 without a warning. Above 250 the conjugate gradient stops when the recurrence residual norm is below 1e-13 times the norm of the right-hand side, or after max(2000, 4 × unknowns) steps, or when pᵀAp is not positive.
water.js
lines 289 to 323
CivilCraft 0.8.1
/** Preconditioned conjugate gradient (Jacobi) on the sparse SPD system; dense Cholesky when small. */
function solveSpd(M, diag, rows, rhs, x0) {
if (M === 0) return [];
if (M <= 250) {
const A = Array.from({ length: M }, () => new Float64Array(M));
for (let i = 0; i < M; i++) { A[i][i] = diag[i]; for (const [j, v] of rows[i]) A[i][j] += v; }
for (let j = 0; j < M; j++) { // Cholesky A = L Lᵀ, in place (lower)
let s = A[j][j]; for (let k = 0; k < j; k++) s -= A[j][k] * A[j][k];
const d = Math.sqrt(Math.max(s, 1e-30)); A[j][j] = d;
for (let i = j + 1; i < M; i++) { let t = A[i][j]; for (let k = 0; k < j; k++) t -= A[i][k] * A[j][k]; A[i][j] = t / d; }
}
const y = new Float64Array(M), x = new Float64Array(M);
for (let i = 0; i < M; i++) { let t = rhs[i]; for (let k = 0; k < i; k++) t -= A[i][k] * y[k]; y[i] = t / A[i][i]; }
for (let i = M - 1; i >= 0; i--) { let t = y[i]; for (let k = i + 1; k < M; k++) t -= A[k][i] * x[k]; x[i] = t / A[i][i]; }
return x;
}
const mul = (v, out) => { for (let i = 0; i < M; i++) { let s = diag[i] * v[i]; for (const [j, a] of rows[i]) s += a * v[j]; out[i] = s; } };
const x = Float64Array.from(x0), r = new Float64Array(M), z = new Float64Array(M), p = new Float64Array(M), Ap = new Float64Array(M);
mul(x, Ap); for (let i = 0; i < M; i++) r[i] = rhs[i] - Ap[i];
for (let i = 0; i < M; i++) { z[i] = r[i] / diag[i]; p[i] = z[i]; }
let rz = 0; for (let i = 0; i < M; i++) rz += r[i] * z[i];
let bn = 0; for (let i = 0; i < M; i++) bn += rhs[i] * rhs[i]; bn = Math.sqrt(bn) || 1;
for (let k = 0; k < Math.max(2000, M * 4); k++) {
let rn = 0; for (let i = 0; i < M; i++) rn += r[i] * r[i];
if (Math.sqrt(rn) < 1e-13 * bn) break;
mul(p, Ap); let pAp = 0; for (let i = 0; i < M; i++) pAp += p[i] * Ap[i];
if (!(pAp > 0)) break;
const al = rz / pAp;
for (let i = 0; i < M; i++) { x[i] += al * p[i]; r[i] -= al * Ap[i]; z[i] = r[i] / diag[i]; }
let rz2 = 0; for (let i = 0; i < M; i++) rz2 += r[i] * z[i];
const be = rz2 / rz; rz = rz2;
for (let i = 0; i < M; i++) p[i] = z[i] + be * p[i];
}
return x;
}
Heads, pressures, velocities and head loss per 1000
64-bit
The numbers shown for each node and pipe, converted back to project units, and the limit checks on them.
pressure = (H − z) × p factor; v = |Q| ÷ (π D² ÷ 4); loss per 1000 = |Hₐ − H_b| ÷ (L ÷ 1000)
Pressure is the difference of total head and elevation, so for a node at high elevation and low pressure it subtracts two close numbers. Flows are divided by the same factor that the demands were multiplied by. Limits are compared with no tolerance.
The largest draw at one node that keeps the lowest pressure in the network at or above the residual pressure.
bisection on the draw: 40 halvings of [lo, hi], each one a full network solve
The upper bound starts at 3000 gpm or 200 L/s and doubles at most 6 times; if the residual pressure is still met at that bound, the halving converges on the bound itself, so the answer is capped at 64 times the starting bound. After 40 halvings the interval is far narrower than the 1e-7 convergence tolerance of each solve, which is what limits the answer.
water.js
lines 325 to 345
CivilCraft 0.8.1
/** Fire flow at one node: the pressures with `flow` drawn there (every other demand as it is), and the flow that brings the
* lowest pressure anywhere to `residual` (bisection on the draw). */
export function fireFlow(water, units, nodeId, { flow, residual, mult = 1 } = {}) {
const W = cleanWater(water, units), d = defaults(units);
const target = residual ?? (units === 'm' ? 140 : 20);
const lowest = (f) => {
const r = solve(W, units, { mult, extra: { [nodeId]: f } });
let m = Infinity, at = null;
for (const [id, n] of r.nodes) if (n.pressure != null && n.pressure < m) { m = n.pressure; at = id; }
return { r, min: m, at };
};
const node = W.nodes.find((n) => n.id === nodeId);
if (!node || node.kind !== 'junction') return { error: 'node' };
const atFlow = flow != null ? lowest(flow) : null;
let lo = 0, hi = units === 'm' ? 200 : 3000, hiR = lowest(hi);
let grow = 0; while (hiR.min >= target && grow++ < 6) { lo = hi; hi *= 2; hiR = lowest(hi); }
const base = lowest(0);
let available = null;
if (base.min >= target) { for (let i = 0; i < 40; i++) { const mid = (lo + hi) / 2; if (lowest(mid).min >= target) lo = mid; else hi = mid; } available = lo; }
return { atFlow: atFlow && { flow, minPressure: atFlow.min, at: atFlow.at, nodePressure: atFlow.r.nodes.get(nodeId).pressure, problems: atFlow.r.problems }, available, residual: target, base: base.min, d };
}
Tank levels over time
64-bit
Steps the network through the day and moves each tank level by the flow into it.
level₊ = level + inflow × Δt ÷ (π d² ÷ 4), clamped to [min level, max level]
Explicit steps: the inflow from one steady state is held constant for the whole time step, so the level history depends on the step length. The number of steps is Math.round(hours × 3600 ÷ step). A level that passes a limit is clamped and reported as a tankfull or tankempty problem.
water.js
lines 347 to 372
CivilCraft 0.8.1
/** Extended period: tank levels and demand pattern over time; hourly (or `step`) hydraulic states. Tanks clamp at their limits and
* the step is reported (no mass is invented: the clamp is a flagged event, a real tank would close its inlet). */
export function extended(water, units, { hours, step } = {}) {
const W = cleanWater(water, units), K = si(units), S = W.settings;
const dt = step || S.step, T = Math.round((hours ?? S.duration) * 3600 / dt);
const levels = {}; for (const n of W.nodes) if (n.kind === 'tank') levels[n.id] = n.level;
const steps = [], problems = [];
for (let s = 0; s <= T; s++) {
const mult = S.pattern.length ? S.pattern[s % S.pattern.length] : 1;
const pumpOn = {}; for (const p of W.pipes) if (p.kind === 'pump' && p.pattern?.length) pumpOn[p.id] = p.pattern[s % p.pattern.length] > 0;
const r = solve(W, units, { mult, levels, pumpOn });
const tanks = {};
for (const n of W.nodes) if (n.kind === 'tank') tanks[n.id] = { level: levels[n.id], inflow: r.nodes.get(n.id).inflow };
steps.push({ t: s * dt, mult, nodes: r.nodes, pipes: r.pipes, tanks, ok: r.ok, problems: r.problems });
if (s === T) break;
for (const n of W.nodes) if (n.kind === 'tank') {
const inflow = (r.nodes.get(n.id).inflow || 0) * K.q; // m³/s into the tank
const areaSI = Math.PI * Math.pow(n.diameter * K.len, 2) / 4;
let lv = levels[n.id] + inflow * dt / areaSI / K.len;
if (lv > n.maxLevel) { problems.push({ code: 'tankfull', id: n.id, t: (s + 1) * dt }); lv = n.maxLevel; }
if (lv < n.minLevel) { problems.push({ code: 'tankempty', id: n.id, t: (s + 1) * dt }); lv = n.minLevel; }
levels[n.id] = lv;
}
}
return { steps, problems, dt };
}
EPANET-format export
64-bit
Writes the network as EPANET-format text for other programs.
written value = round(v × 10⁴) ÷ 10⁴
Every number written (coordinates, elevations, demands, lengths, diameters, settings, curve points, pattern values) is rounded to 4 decimals by Math.round(v × 1e4) ÷ 1e4 before it is turned into text. The model itself is not changed, but a file read back carries the rounded values.
water.js
lines 378 to 402
CivilCraft 0.8.1
export function writeInp(water, units) {
const W = cleanWater(water, units), M = units === 'm', f = (v, d = 4) => (v == null ? '' : String(Math.round(v * 10 ** d) / 10 ** d));
const L = ['[TITLE]', 'CadCraft water mains', '', '[JUNCTIONS]', ';ID Elevation Demand'];
for (const n of W.nodes) if (n.kind === 'junction') L.push(`${n.id}\t${f(n.z ?? 0)}\t${f(n.demand)}`);
L.push('', '[RESERVOIRS]', ';ID Head');
for (const n of W.nodes) if (n.kind === 'reservoir') L.push(`${n.id}\t${f(n.head ?? n.z ?? 0)}`);
L.push('', '[TANKS]', ';ID Elevation InitLevel MinLevel MaxLevel Diameter MinVol');
for (const n of W.nodes) if (n.kind === 'tank') L.push(`${n.id}\t${f(n.z ?? 0)}\t${f(n.level)}\t${f(n.minLevel)}\t${f(n.maxLevel)}\t${f(n.diameter)}\t0`);
L.push('', '[PIPES]', ';ID Node1 Node2 Length Diameter Roughness MinorLoss Status');
const pipeLen = (p) => p.length ?? Math.hypot(W.nodes.find((n) => n.id === p.a).x - W.nodes.find((n) => n.id === p.b).x, W.nodes.find((n) => n.id === p.a).y - W.nodes.find((n) => n.id === p.b).y);
for (const p of W.pipes) if (p.kind === 'pipe') L.push(`${p.id}\t${p.a}\t${p.b}\t${f(pipeLen(p))}\t${f(p.d)}\t${f(p.c ?? W.settings.c)}\t${f(p.minor)}\t${p.status === 'closed' ? 'Closed' : p.check ? 'CV' : 'Open'}`);
L.push('', '[PUMPS]', ';ID Node1 Node2 Parameters');
const curves = [];
for (const p of W.pipes) if (p.kind === 'pump') { curves.push(p); L.push(`${p.id}\t${p.a}\t${p.b}\tHEAD C_${p.id}`); }
L.push('', '[VALVES]', ';ID Node1 Node2 Diameter Type Setting MinorLoss');
for (const p of W.pipes) if (p.kind === 'prv') L.push(`${p.id}\t${p.a}\t${p.b}\t${f(p.d)}\tPRV\t${f(p.setting)}\t${f(p.minor)}`);
L.push('', '[CURVES]', ';ID Flow Head');
for (const p of curves) for (const [q, h] of p.curve) L.push(`C_${p.id}\t${f(q)}\t${f(h)}`);
if (W.settings.pattern.length) { L.push('', '[PATTERNS]'); for (let i = 0; i < W.settings.pattern.length; i += 6) L.push('1\t' + W.settings.pattern.slice(i, i + 6).map((v) => f(v)).join('\t')); }
L.push('', '[TIMES]', `Duration\t${W.settings.duration}:00`, `Hydraulic Timestep\t${Math.round(W.settings.step / 3600 * 100) / 100 >= 1 ? Math.round(W.settings.step / 3600) + ':00' : '0:' + String(Math.round(W.settings.step / 60)).padStart(2, '0')}`);
L.push('', '[OPTIONS]', `Units\t${M ? 'LPS' : 'GPM'}`, 'Headloss\tH-W', M ? 'Pattern\t1' : 'Pattern\t1', '', '[COORDINATES]', ';Node X Y');
for (const n of W.nodes) L.push(`${n.id}\t${f(n.x)}\t${f(n.y)}`);
L.push('', '[END]', '');
return L.join('\n');
}
In detail
Method
Global gradient method (the file header cites Todini and Pilati 1988) with Hazen-Williams pipe loss, minor loss K·v²/2g and power-law pump curves. Each pass solves one symmetric positive definite system: dense Cholesky for 250 unknowns or fewer, otherwise conjugate gradient with a Jacobi preconditioner. Fire flow is a bisection on the draw; extended period is one steady state per time step with explicit tank-level updates.
Number types
64-bit numbers and Float64Array. No 32-bit arrays.
Tolerances
Outer iteration: Σ|ΔQ| ÷ Σ|Q| below 1e-7 with no status change, at most 300 passes. Conjugate gradient: residual norm below 1e-13 × the norm of the right-hand side, at most max(2000, 4 × unknowns) steps. Flow floor in the loss terms 1e-7 m³/s. Status switches use 1e-9 m³/s for reverse flow and 1e-6 m or 1e-9 m for head comparisons, all absolute. A Cholesky pivot is floored at √1e-30.
Limits
Hazen-Williams is an empirical formula with a coefficient of 10.67 and exponents 1.852 and 4.871 as written. GPM = 6.30902e-5 m³/s and 0.433527 psi per foot are rounded constants. A closed link stays in the matrix with a conductance of 1e-8. The EPANET-format export rounds every number to 4 decimals.
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.