Resources Numerical accuracy Calculation policy
Runs rain over a grid of the graded ground and reports where water ponds, how deep it gets, what overtops, and the volumes that fall, soak in, drain and leave.
The simulation state (depths, face flows, volumes) is 64-bit and standard gravity 9.80665 is used. The results are set by the method: an explicit time-stepping scheme on a grid. The per-cell maximum depth, maximum speed and first-wet time are kept in 32-bit arrays, and the deepest water shown in the results is the maximum of the 32-bit depth array.
Open item. Keep the per-cell maximum depth (and the reported deepest value) in 64-bit, or track the overall maximum in a 64-bit variable, instead of reading it back from the Float32Array.
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.
Tests bin/test-civilcraft-rain.mjs
The constants the simulation uses.
g = 9.80665 m/s²; 1 ft = 0.3048 m
Both are exact defined values. The foot here is the international foot, while rail.js and the IFC export convert a feet project with the US survey foot.
export const G = 9.80665;
const FT = 0.3048;Samples the ground surface into square cells in meters.
n_x = round(W ÷ c), n_y = round(H ÷ c); c = max(W ÷ n_x, H ÷ n_y); z (m) = elevation × meters per unit
Cell centers come from whole-number indexes. The cell is stretched to the larger of the two fits, so the grid can reach slightly past the bounds on one side; cells with no ground are stored as NaN. Elevations are multiplied by 0.3048 for feet projects and kept in a Float64Array.
export function demFromSurface(surface, { cell, units = 'ft', bounds = null, maxCells = 4e6, mask = null } = {}) {
if (!surface) return null;
const b = bounds || surface.stats().bounds;
const W = b.maxX - b.minX, H = b.maxY - b.minY;
if (!(W > 0 && H > 0)) return null;
let c = cell > 0 ? cell : Math.sqrt((W * H) / 40000);
while ((W / c) * (H / c) > maxCells) c *= 1.25;
const unit = units === 'm' ? 1 : 1 / FT, f = units === 'm' ? 1 : FT; // unit = project length per metre; f = metres per project length
const nx = Math.max(1, Math.round(W / c)), ny = Math.max(1, Math.round(H / c));
c = Math.max(W / nx, H / ny);
const z = new Float64Array(nx * ny);
for (let j = 0; j < ny; j++) for (let i = 0; i < nx; i++) {
const x = b.minX + (i + 0.5) * c, y = b.minY + (j + 0.5) * c;
const v = surface.elevationAt(x, y);
z[j * nx + i] = v === null || v === undefined ? NaN : v * f;
}
const dem = { nx, ny, cell: c, cellM: c * f, unit, x0: b.minX, y0: b.minY, z, block: new Uint8Array(nx * ny), n: null, loss: null };
if (mask) for (let k = 0; k < z.length; k++) if (!mask(dem.x0 + ((k % nx) + 0.5) * c, dem.y0 + (Math.floor(k / nx) + 0.5) * c)) z[k] = NaN;
return dem;
}Finds each depression, its area, the volume it holds before it spills, its greatest depth and its sill level.
filled = max(ground, lowest way out); depth = filled − ground; volume = Σ depth × cell²
Priority-Flood on 64-bit elevations with a binary heap; equal keys are compared exactly. A cell belongs to a pond when its fill depth is over the minimum (0.005 m by default). Volume is a running sum of depths times the cell area.
export function findDepressions(dem, { minDepth = 0.005 } = {}) {
const { nx, ny, z, block } = dem, N = nx * ny;
const F = new Float64Array(N).fill(Infinity), done = new Uint8Array(N), heap = new Heap();
const usable = (k) => !Number.isNaN(z[k]) && !block[k];
const nb = [[-1, 0], [1, 0], [0, -1], [0, 1]];
for (let j = 0; j < ny; j++) for (let i = 0; i < nx; i++) {
const k = j * nx + i;
if (!usable(k)) continue;
let edge = i === 0 || j === 0 || i === nx - 1 || j === ny - 1;
if (!edge) for (const [di, dj] of nb) if (Number.isNaN(z[k + dj * nx + di])) { edge = true; break; } // next to no ground = leaves the map
if (edge) { F[k] = z[k]; done[k] = 1; heap.push(z[k], k); }
}
while (heap.size) {
const k = heap.pop(), i = k % nx, j = (k - i) / nx;
for (const [di, dj] of nb) {
const a = i + di, b = j + dj;
if (a < 0 || b < 0 || a >= nx || b >= ny) continue;
const m = b * nx + a;
if (done[m] || !usable(m)) continue;
done[m] = 1;
F[m] = Math.max(z[m], F[k]);
heap.push(F[m], m);
}
}
const label = new Int32Array(N).fill(-1), ponds = [];
const stack = [];
for (let k = 0; k < N; k++) {
if (label[k] >= 0 || !done[k] || F[k] - z[k] <= minDepth) continue;
const id = ponds.length, p = { id, cells: 0, area: 0, volume: 0, maxDepth: 0, sill: -Infinity, x: 0, y: 0, low: z[k] };
stack.push(k); label[k] = id;
while (stack.length) {
const q = stack.pop(), i = q % nx, j = (q - i) / nx, d = F[q] - z[q];
p.cells++; p.volume += d; if (d > p.maxDepth) p.maxDepth = d; if (F[q] > p.sill) p.sill = F[q];
if (z[q] < p.low) { p.low = z[q]; p.x = dem.x0 + (i + 0.5) * dem.cell; p.y = dem.y0 + (j + 0.5) * dem.cell; }
for (const [di, dj] of nb) {
const a = i + di, b = j + dj;
if (a < 0 || b < 0 || a >= nx || b >= ny) continue;
const m = b * nx + a;
if (label[m] >= 0 || !done[m] || F[m] - z[m] <= minDepth) continue;
label[m] = id; stack.push(m);
}
}
p.area = p.cells * dem.cellM * dem.cellM; p.volume *= dem.cellM * dem.cellM;
ponds.push(p);
}
return { filled: F, label, ponds };
}The total depth of rain over the storm's duration.
depth (mm) = Σ i × Δt ÷ 60
A sum of intensity times block length in minutes; 60 is exact.
export function hyetoDepth(rain, duration) {
if (typeof rain === 'number') return (rain * duration) / 60;
let mm = 0;
for (let k = 0; k < rain.length; k++) { const a = rain[k].t, b = k + 1 < rain.length ? Math.min(rain[k + 1].t, duration) : duration; if (b > a) mm += (rain[k].i * (b - a)) / 60; }
return mm;
}Builds a storm from an intensity-duration curve with the heaviest block in the middle.
depth_k = idf(k·Δt) × k·Δt ÷ 60 − depth_(k−1); i_k = depth_k × 60 ÷ Δt
Each block is the difference of two cumulative depths, so for long durations it is a small difference of larger numbers. The number of blocks is the duration over the step, rounded.
/** Alternating blocks: the biggest block in the middle, the next to its right, then left, and so on. */
export function alternatingBlock(idf, duration, step = 5) {
const n = Math.max(1, Math.round(duration / step)), depths = [];
let prev = 0;
for (let k = 1; k <= n; k++) { const d = (idf(k * step) * k * step) / 60; depths.push(d - prev); prev = d; }
const out = new Array(n); let lo = Math.floor((n - 1) / 2), hi = lo + 1, right = true;
out[lo] = depths[0]; lo--;
for (let k = 1; k < n; k++) { if (right && hi < n) out[hi++] = depths[k]; else if (lo >= 0) out[lo--] = depths[k]; else out[hi++] = depths[k]; right = !right; }
return out.map((d, k) => ({ t: k * step, i: (d * 60) / step }));
}The length of each simulation step.
dt = min(dt_max, α × Δx ÷ √(g × max(h_max, 0.01)))
The step follows the deepest water of the step before, with a floor of 0.01 m on that depth and a cap of 5 s by default. Time advances by adding dt each step, so the clock is a running sum.
let dt = Math.min(o.dtMax, (o.alpha * dx) / Math.sqrt(G * Math.max(this.hmax, 0.01)));
const tEndS = o.tEnd * 60;
if (this.t + dt > tEndS) dt = tEndS - this.t;
if (dt <= 0) { this.done = true; return; }The unit flow between two neighboring cells for one time step.
q = (q − g·h_f·dt·Δη ÷ Δx) ÷ (1 + g·dt·n²·|q| ÷ h_f^(7/3)); |q| ≤ h_f·√(g·h_f)
One 64-bit expression per face with Math.pow for the 7/3 power. The flow depth is the higher water level minus the higher bed; faces with 1e-5 m or less carry nothing, which also keeps the power from dividing by zero. The flow is capped at a Froude number of 1.
const ea = z[a] + h[a], eb = z[b] + h[b], hf = Math.max(ea, eb) - Math.max(z[a], z[b]);
if (hf <= hmin) { qx[f] = 0; continue; }
const n = (nOf(a) + nOf(b)) / 2, q0 = qx[f];
let q = (q0 - G * hf * dt * ((eb - ea) / dx)) / (1 + (G * dt * n * n * Math.abs(q0)) / Math.pow(hf, 7 / 3));
const lim = hf * Math.sqrt(G * hf); if (q > lim) q = lim; else if (q < -lim) q = -lim;
qx[f] = q;The unit flow out of the grid at an open edge or beside a hole in the ground.
q = h^(5/3) × √S ÷ n, S = max(edge slope, bed slope to the cell behind); q ≤ h·√(g·h)
Manning's formula for a wide sheet of water in SI units, with a default edge slope of 0.001 where the bed gives none.
let S = o.edgeSlope;
if (inner >= 0 && wet(inner)) S = Math.max(S, (z[c] - z[inner]) / dx);
const q = Math.pow(h[c], 5 / 3) * Math.sqrt(S) / nOf(c), lim = h[c] * Math.sqrt(G * h[c]);
qx[f] = (A ? 1 : -1) * Math.min(q, lim);Scales the outflows of a cell so it cannot give away more water than it has.
when Σ out × dt ÷ Δx > available: out ← out × available ÷ (Σ out × dt ÷ Δx)
The scale factor is a quotient of two 64-bit sums. A face shared with a neighbor can be scaled from either side; the update that follows uses the scaled values, so the mass leaving one cell is the mass entering the next.
const avail = Math.max(0, h[k] + add); // the water this cell may give away this step
if (out * k2 > avail && out > 0) {
const s = avail / (out * k2);
if (qx[fW] < 0) qx[fW] *= s; if (qx[fE] > 0) qx[fE] *= s; if (qy[fS] < 0) qy[fS] *= s; if (qy[fN] > 0) qy[fN] *= s;
}Updates each cell's depth and records its greatest depth, greatest speed and the time it first got wet.
h ← h + rain − losses + (dt ÷ Δx) × (q_W − q_E + q_S − q_N); speed = |q| ÷ h
The depth itself is 64-bit, but the three records are Float32Array: each new maximum is compared with the stored 32-bit value and then rounded to 32-bit when stored. The deepest water shown in the results is the maximum of this 32-bit depth array, so it carries about 7 digits. A depth below zero is set to 0, and flagged when below −1e-12.
const inflow = qx[fW] - qx[fE] + qy[fS] - qy[fN];
let hn = h[k] + tmp[k] + k2 * inflow;
if (hn < 0) { if (hn < -1e-12) this.negativeSeen = true; hn = 0; }
h[k] = hn;
if (hn > hmax) hmax = hn;
if (hn > this.maxDepth[k]) this.maxDepth[k] = hn;
if (hn > hmin) {
if (this.firstWet[k] < 0) this.firstWet[k] = this.t + dt;
const u = (qx[fW] + qx[fE]) / 2, v = (qy[fS] + qy[fN]) / 2, sp = Math.hypot(u, v) / hn;
if (sp > this.maxSpeed[k]) this.maxSpeed[k] = sp;Sets up the arrays that hold the per-cell maximum depth, maximum speed and first-wet time.
All three are created as Float32Array. The saved frames and the end depth are also 32-bit copies of the 64-bit depth; those two are used to draw the map.
this.maxDepth = new Float32Array(N); this.maxSpeed = new Float32Array(N); this.firstWet = new Float32Array(N).fill(-1);Takes water out of a drain's cell up to its capacity each step.
take = min(capacity × dt, depth × cell area)
The capacity was converted from liters per second to cubic meters per second by the exact 1000 when the drain was set up. The volume taken and offered accumulate as 64-bit running sums.
const want = d.cap * dt, have = h[d.cell] * this.A, take = Math.min(want, have);
d.offered += Math.max(want, 0); d.taken += take; this.vol.drained += take; h[d.cell] -= take / this.A;
if (h[d.cell] < 0) h[d.cell] = 0;Adds each step's rain, losses and outflow to the running totals.
out += Σ q × Δx × dt; rain += Σ depth of rain × cell area
The totals are plain 64-bit running sums over all steps, with no compensation.
const outStep = outV * dx * dt; // m³ this step (q is m²/s along a face of length dx)
this.vol.out += outStep; this.vol.rain += rainV * this.A; this.vol.loss += lossV * this.A;The stored volume, the mass balance error, the pond volumes at the end and the peak outflow.
error = (rain − loss − out − drained − stored) ÷ rain
The error is a difference of large totals divided by the rain volume; a result with an error above 2 percent is flagged. The peak outflow is the largest recorded value of the hydrograph, which is sampled about 400 times over the run, not at every step.
result() {
const { dem, h } = this, N = dem.nx * dem.ny;
let stored = 0; for (let m = 0; m < this.live.length; m++) stored += h[this.live[m]];
stored *= this.A;
const v = this.vol, error = v.rain > 0 ? (v.rain - v.loss - v.out - v.drained - stored) / v.rain : 0;
const F = this._filled || (this._filled = findDepressions(dem, { minDepth: this.o.minPond }).filled);
const ponds = this.ponds.map((p) => {
let held = 0; for (const k of p.cellList) held += Math.min(h[k], Math.max(0, F[k] - dem.z[k]) + 0) ;
return { id: p.id, cells: p.cells, area: p.area, volume: p.volume, maxDepth: p.maxDepth, sill: p.sill, x: p.x, y: p.y, peakVolume: p.peakVolume, endVolume: held * this.A, overtopped: p.overtopped, overtopAt: p.overtopAt, spilled: this.spillOut[p.id] };
});
let peak = { q: 0, t: 0 }; for (const r of this.hydro) if (r.out > peak.q) peak = { q: r.out, t: r.t };
const problems = [];
if (this.negativeSeen) problems.push({ code: 'negative' });
if (Math.abs(error) > 0.02) problems.push({ code: 'balance', error });
for (const p of ponds) if (p.overtopped) problems.push({ code: 'overtop', pond: p.id, x: p.x, y: p.y, at: p.overtopAt, spilled: p.spilled });
for (const d of this.drains) if (d.offered > 0 && d.taken < d.offered * 0.5 && d.taken < this.vol.rain * 0.0 + 1e9) { /* a drain that is not asked for its capacity is fine */ }
return {
t: this.t, steps: this.steps, cells: this.activeCells, volumes: { rain: v.rain, loss: v.loss, out: v.out, drained: v.drained, stored, error },
hydrograph: this.hydro, maxDepth: this.maxDepth, maxSpeed: this.maxSpeed, firstWet: this.firstWet, endDepth: Float32Array.from(h), frames: this.frames,
ponds, drains: this.drains.map((d) => ({ name: d.name, taken: d.taken, offered: d.offered, bypass: Math.max(0, d.offered - d.taken) })), gauges: this.gauges, peak, problems, dem: { nx: dem.nx, ny: dem.ny, cell: dem.cell, cellM: dem.cellM, unit: dem.unit, x0: dem.x0, y0: dem.y0 },
};
}One step stores or exchanges 32-bit values that reach a number a person reads. It is listed here with its open item, and the step is either moved to 64-bit or the app will warn where the number is read.