Routes a flood hydrograph down the cross-sections of a river model and reports water levels, flows, their peaks and a volume balance over time.
How the precision is kept
The state (water level and flow at each section) is held in Float64Array. Accuracy is set by the scheme: the time step and section spacing, area and conveyance read from a 201-level table by linear interpolation, a forward-difference Jacobian, and Newton iteration stopped at a relative change of 1e-8. A step that does not converge is halved up to 4 times. The result carries its own volume balance error.
Open item. Write the exponent as 5 ÷ 3 instead of 1.6667 where conveyance is extended below the lowest table level.
Calculations under this policy
6 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.
Flow area and conveyance of one cross-section as functions of water level, read from a table.
levels hᵢ = h₀ + (hi − h₀)·(i ÷ 200)²; between levels: linear interpolation; below h₀: A·f, K·f^1.6667 with f = depth ÷ (h₀ − bed)
The table has 201 levels spaced as the square of the index, and values between levels are linear interpolations, so area and conveyance are approximations of the exact section values. Below the lowest level conveyance is scaled by the exponent 1.6667, a rounded 5/3. Above the top level the last interval is extended in a straight line.
river-unsteady.js
lines 33 to 49
CivilCraft 0.8.1
/** A(h) and K(h) of one section as a table of LEVELS + 1 levels, closer together near the bed (levels grow as the square of the index). */
function table(geo, units) {
const lo = geo.minZ, depth = Math.max(geo.maxZ - geo.minZ, 1e-3), h0 = lo + 0.005 * depth, hi = geo.maxZ + depth;
const h = [], A = [], K = [];
for (let i = 0; i <= LEVELS; i++) {
const w = h0 + (hi - h0) * (i / LEVELS) ** 2, r = hydraulics(geo, w, units);
h.push(w); A.push(r.A); K.push(r.K);
}
const at = (arr, w) => {
if (w <= h0) { const f = Math.max(w - lo, 1e-9) / (h0 - lo); return arr[0] * (arr === A ? f : Math.pow(f, 1.6667)); }
if (w >= hi) { const a = arr[LEVELS], b = arr[LEVELS - 1]; return a + ((w - hi) / (h[LEVELS] - h[LEVELS - 1])) * (a - b); }
const i = Math.min(LEVELS - 1, Math.floor(LEVELS * Math.sqrt((w - h0) / (hi - h0))));
const j = w < h[i] ? i - 1 : w > h[i + 1] ? i + 1 : i, k = Math.max(0, Math.min(LEVELS - 1, j));
return arr[k] + ((w - h[k]) / (h[k + 1] - h[k])) * (arr[k + 1] - arr[k]);
};
return { lo, h0, A: (w) => at(A, w), K: (w) => at(K, w) };
}
Banded linear solve
64-bit
Solves the Newton step for all sections at once.
L·U·x = b with row pivoting inside the band
Gaussian elimination in 64-bit with partial pivoting over the two rows below the diagonal; the solution is a Float64Array. A pivot of magnitude 1e-300 or less returns no solution. No iterative refinement is done; the Newton loop that calls it repeats until the step is small.
river-unsteady.js
lines 58 to 83
CivilCraft 0.8.1
/** Banded LU solve (partial pivoting; kl = 2, ku = 3, fill up to 5 to the right): rows as short arrays on an absolute column window. */
function bandSolve(rows, c0, rhs) {
const n = rhs.length, W = rows[0].length, KL = 2, REACH = 6;
const get = (r, c) => { const j = c - c0[r]; return j >= 0 && j < W ? rows[r][j] : 0; };
const sub = (r, c, v) => { const j = c - c0[r]; if (j < 0 || j >= W) { if (v === 0) return; throw new Error('band overflow'); /* i18n-ignore: internal */ } rows[r][j] -= v; };
for (let k = 0; k < n; k++) {
let p = k, best = Math.abs(get(k, k));
for (let i = k + 1; i <= Math.min(n - 1, k + KL); i++) { const v = Math.abs(get(i, k)); if (v > best) { best = v; p = i; } }
if (!(best > 1e-300)) return null;
if (p !== k) { [rows[k], rows[p]] = [rows[p], rows[k]]; [c0[k], c0[p]] = [c0[p], c0[k]]; [rhs[k], rhs[p]] = [rhs[p], rhs[k]]; }
const piv = get(k, k);
for (let i = k + 1; i <= Math.min(n - 1, k + KL); i++) {
const f = get(i, k) / piv;
if (f === 0) continue;
for (let c = k; c <= Math.min(n - 1, k + REACH); c++) sub(i, c, f * get(k, c));
rhs[i] -= f * rhs[k];
}
}
const x = new Float64Array(n);
for (let k = n - 1; k >= 0; k--) {
let s = rhs[k];
for (let c = k + 1; c <= Math.min(n - 1, k + REACH); c++) s -= get(k, c) * x[c];
x[k] = s / get(k, k);
}
return x;
}
Continuity and momentum residuals of one interval
64-bit
The two equations that must be zero between a pair of neighboring sections for one time step.
64-bit sums of differences between new and old values. When the flow changes little in a step, the time terms subtract nearly equal sums, which is inherent to the scheme. Division by area is guarded (zero area gives zero); conveyance is floored at 1e-12 before it is squared.
Finds the water levels and flows at the end of a time step.
J·Δu = −F; J by forward differences with steps 1e-6·Href and 1e-6·Qref; stop when max(|Δh| ÷ Href, |ΔQ| ÷ Qref) < 1e-8
The Jacobian is a forward difference, so it is approximate; the residual itself is evaluated exactly, and the loop stops on the size of the update relative to the reference depth and flow (1e-8), within 20 iterations. Level updates are capped at half the reference depth and levels are held at least 1e-4 × reference depth above the bed.
river-unsteady.js
lines 137 to 178
CivilCraft 0.8.1
const attempt = (tNew, dt) => {
const ho = h, Qo = Q;
let hn = Float64Array.from(ho), Qn = Float64Array.from(Qo);
const n = 2 * N;
for (let it = 1; it <= 20; it++) {
const F = new Float64Array(n);
const rows = [], c0 = [];
const addRow = (r) => { rows[r] = new Array(12).fill(0); c0[r] = Math.max(0, r - 3); };
for (let r = 0; r < n; r++) addRow(r);
const put = (r, c, v) => { rows[r][c - c0[r]] = v; };
// upstream: Q0 = inflow(t)
F[0] = Qn[0] - linear(inflow, tNew, 'q'); put(0, 1, 1);
for (let i = 0; i + 1 < N; i++) {
const base = [hn[i], Qn[i], hn[i + 1], Qn[i + 1]], r0 = interval(i, ...base, ho[i], Qo[i], ho[i + 1], Qo[i + 1], dt);
F[2 * i + 1] = r0[0]; F[2 * i + 2] = r0[1];
for (let v = 0; v < 4; v++) {
const b = base.slice(); const e = v % 2 === 0 ? dhEps : dqEps; b[v] += e;
const r1 = interval(i, ...b, ho[i], Qo[i], ho[i + 1], Qo[i + 1], dt), col = 2 * i + v;
put(2 * i + 1, col, (r1[0] - r0[0]) / e); put(2 * i + 2, col, (r1[1] - r0[1]) / e);
}
}
const L = N - 1;
if (down.type === 'normal') {
const k = Kf(L, hn[L]), sl = Math.sqrt(down.slope), e = dhEps;
F[n - 1] = Qn[L] - k * sl; put(n - 1, 2 * L, -((Kf(L, hn[L] + e) - k) / e) * sl); put(n - 1, 2 * L + 1, 1);
} else { F[n - 1] = hn[L] - linear(down.series, tNew, 'h'); put(n - 1, 2 * L, 1); }
const rhs = Float64Array.from(F, (v) => -v);
let du;
try { du = bandSolve(rows, c0, rhs); } catch { du = null; }
if (!du) return null;
let worst = 0;
for (let i = 0; i < N; i++) {
let dh = du[2 * i], dq = du[2 * i + 1];
const cap = 0.5 * Href; if (Math.abs(dh) > cap) dh = Math.sign(dh) * cap;
hn[i] = Math.max(tabs[i].lo + 1e-4 * Href, hn[i] + dh); Qn[i] += dq;
worst = Math.max(worst, Math.abs(dh) / Href, Math.abs(dq) / Qref);
}
newtonMax = Math.max(newtonMax, it);
if (worst < 1e-8) return { h: hn, Q: Qn };
}
return null;
};
Time stepping, peaks and Froude check
64-bit
Advances the run step by step, records the peak level and flow at each section and flags supercritical flow.
V_in += Δt·(θ·Q_new + (1 − θ)·Q_old); T ≈ (A(h + 0.001·Href) − A(h)) ÷ (0.001·Href); flag when |Q| ÷ A > √(g·A ÷ T)
Inflow and outflow volumes are running 64-bit sums, one term per step, added in time order with no compensation. Peaks are compared and stored in 64-bit. The top width for the Froude check is a finite difference of the interpolated area table. A failed step is retried with half the time step, up to 4 times.
river-unsteady.js
lines 179 to 196
CivilCraft 0.8.1
while (t < dur - 1e-9) {
let dt = Math.min(dt0, dur - t), res = attempt(t + dt, dt), tries = 0;
while (!res && tries < 4) { dt /= 2; tries++; halved++; res = attempt(t + dt, dt); }
if (!res) { flag('noconv'); break; }
// volumes at the two ends, weighted in time like the scheme itself
const L = N - 1;
mass.inflow += dt * (theta * res.Q[0] + (1 - theta) * Q[0]);
mass.outflow += dt * (theta * res.Q[L] + (1 - theta) * Q[L]);
h = res.h; Q = res.Q; t += dt; steps++;
for (let i = 0; i < N; i++) {
if (h[i] > peaks[i].wse) { peaks[i].wse = h[i]; peaks[i].tWse = t; }
if (Q[i] > peaks[i].q) { peaks[i].q = Q[i]; peaks[i].tQ = t; }
const a = Af(i, h[i]), top = geos[i].maxZ;
if (h[i] > top) flag('walls');
if (a > 0) { const T = (Af(i, h[i] + 1e-3 * Href) - a) / (1e-3 * Href); if (T > 0 && Math.abs(Q[i]) / a > Math.sqrt((g * a) / T)) flag('froude'); }
}
if (t >= nextFrame - 1e-9 || t >= dur - 1e-9) { frames.push({ t, h: Array.from(h), q: Array.from(Q) }); nextFrame += every; }
}
Volume balance error
64-bit
The share of water not accounted for between inflow, outflow and the change of storage in the reach.
Storage is a trapezoid sum of table areas along the reach. The error subtracts volumes of similar size, so it shows the scheme's balance down to the rounding of those sums; the divisor is floored at 1e-12.
Preissmann four-point implicit scheme for the Saint-Venant equations; the file cites Preissmann 1961 and Cunge, Holly and Verwey 1980, chapter 3. Time weighting θ defaults to 0.6. Newton's method per step with a banded LU solve with row pivoting.
Number types
64-bit numbers and Float64Array.
Tolerances
Newton stops when the largest change is below 1e-8 of the reference depth and reference flow, with at most 20 iterations. Jacobian steps are 1e-6 × reference depth and 1e-6 × reference flow. A pivot at or below 1e-300 fails the solve. Level changes per iteration are capped at half the reference depth; levels are kept at least 1e-4 × reference depth above the bed. Conveyance is floored at 1e-12. Time comparisons use 1e-9 s.
Limits
Subcritical flow only, no structures inside the reach, no lateral inflow, one reach, energy coefficient 1 in the momentum term. Below the lowest table level conveyance is scaled with the exponent 1.6667, a rounded 5/3. Gravity in feet is the rounded 32.174.
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.