01 Quick start

FluxCell solves the two-dimensional shallow-water equations on a uniform grid with WebGPU compute shaders. It is plain ES modules with no dependencies and no build step: serve the src/ folder and import it.

import { compileKernels, createSolver } from './fluxcell/src/index.mjs';

const adapter = await navigator.gpu.requestAdapter();
const device = await adapter.requestDevice();
const kernels = await compileKernels(device, 'muscl-fused'); // compile once

// 200 m x 100 m flat channel, 1 m cells: a dam at x = 50 m, 2 m of water
// behind it and 0.5 m beyond.
const width = 200, height = 100, n = width * height;
const bed = new Float32Array(n);           // bed elevation (m)
const state = new Float32Array(4 * n);     // per cell [h, hu, hv, 0]
for (let i = 0; i < n; i++) state[4 * i] = i % width < 50 ? 2 : 0.5;

const solver = await createSolver(device,
    { width, height, dx: 1, dy: 1,
      boundaryX: 'outflow', boundaryY: 'wall' },
    { state, bed, manning: new Float32Array(n).fill(0.03) },
    kernels);
try {
    await solver.advance(30);              // 30 s of simulated time
    const result = await solver.readState();
    console.log('centre depth', result[4 * (50 * width + 100)], 'm');
} finally {
    solver.destroy();                      // the caller owns the device
}

Rows run south to north and columns west to east; cell i = y * width + x. Depths are in metres, discharges per unit width in m²/s. With a bundler or Node, the package name is fluxcell.

02 Configuration and inputs

createSolver(device, config, initial, kernels) validates everything and rejects inconsistent input (negative or nonfinite depths, momentum in dry cells, water in solid cells) instead of correcting it.

Configuration

Initial data and forcing

// A river reach: 40 m³/s entering on the west edge (ramped over 10 min), the east edge held at 2.5 m.
hydrographs: [
    { edge: 'west', kind: 'discharge', series: [[0, 0], [600, 40]] },
    { edge: 'east', kind: 'stage', series: [[0, 2.5]] },
]

03 Running and control

Host fields: time, sourceVolume, counters (flux fallbacks, order reductions, retries) and allocatedBytes. The double-precision CPU solvers MusclReferenceSolver and ReferenceSolver take the same configuration and run in Node or the browser, for testing and for checking GPU results.

Full reference: docs/API.md.

04 Equations

The depth-averaged shallow-water (Saint-Venant) equations for depth h and discharges hu, hv over a bed z:

and likewise for hv. R and S are rain and sources, the second term on the right is Manning friction, τ the kinematic wind stress and f the Coriolis parameter. A tracer of concentration c obeys

with diffusivity D (plus the eddy diffusivity of the closure, if on) and settling velocity w_s.

05 Numerical scheme

A conservative finite-volume method: each cell holds averages of h, hu, hv, and each step exchanges fluxes across the four faces of every cell, so water is conserved to rounding.

1. Reconstruction (MUSCL). Slopes of surface elevation and velocity in each cell, limited by the monotonized-central limiter, give second-order values at the faces without new extrema. Reconstructing the surface (h + z) and not the depth keeps a lake at rest exactly at rest.
2. Hancock predictor. The face values are advanced by half a step with the local equations, making the scheme second order in time in a single stage.
3. Hydrostatic reconstruction (Audusse et al., 2004) at each face keeps depths nonnegative and balances the bed-slope term against the pressure flux (well-balanced).
4. HLLC Riemann solver (Toro et al., 1994). The flux through each face is the solution of the local Riemann problem with the contact wave restored, which keeps shear layers and tracer fronts sharp. Wet/dry faces and degenerate states fall back to HLLE.
5. Update and sources. The conservative update, then rain and sources, wind stress, the exact Coriolis rotation, and Manning friction applied implicitly (unconditionally stable, never reverses the flow).
6. Wetting and drying. Thin layers below 0.1 mm reconstruct at first order and have their velocity regularized, so fronts advance over dry ground without negative depths or spurious speeds.
7. Timestep. CFL 0.4 from the largest wave speed |u| + √(gh) and the boundary ghost states, computed on the GPU at every step.

A first-order Rusanov scheme is kept as a control.

06 GPU implementation

Three kernels share one interface: muscl-fused (the production path), muscl-hancock (the same method in five auditable passes) and rusanov (first-order control). The fused kernel loads a 20 × 20 tile with a two-cell halo into workgroup memory and performs reconstruction, fluxes, the update and the next step's wave speeds in one dispatch; faces on tile seams are computed identically by both neighbours, so mass still cancels exactly.

• The timestep is chosen on the GPU; the host reads back once per batch of up to 256 steps.
• Depth carries a compensation term (Kahan summation) and the bed is shifted by a double-precision datum, so f32 holds lakes at rest exactly and keeps precision over thousands of metres of relief.
• An invalid step anywhere rolls the batch back to its checkpoint and replays it one step at a time, retrying at half the step.
• Results do not depend on batch size: the same run in batches of 1, 16 or 256 steps is bitwise identical.

07 Verification

Every GPU kernel is compared, step by step, with a double-precision CPU implementation of the same method, and both are compared with exact or measured results.

Details: docs/VERIFICATION.md; the method in full: docs/METHOD.md. Run the suites yourself at tools/validate/ (GPU) and with npm test (CPU).

08 Requirements and deployment

• Browser: WebGPU (Chrome and Edge 113+, recent Chrome on Android, Safari 26+, Firefox 141+ on Windows). The page must be served over HTTPS or from localhost.
• Dependencies: none. The solver and these demos are native ES modules and WGSL; no build step or bundler.
• Hosting: any static server. Upload src/, examples/viewer/ and data/dem/ together (the demos load them by relative path) and serve .mjs as text/javascript.
• Tests: npm test runs the CPU suite in Node 20+; npm run serve serves the repository at http://localhost:8080/ for the GPU suites and these demos.
• Memory: about 100 bytes of GPU memory per cell for the fused kernel with a tracer; a 1,000 × 1,000 grid fits easily on integrated GPUs.

09 Limits

• Shallow-water physics: no vertical structure or frequency dispersion (short waves over deep water, the steep fronts of solitary waves, the 3D horseshoe vortex at a pier).
• Single precision on the GPU (with the safeguards above).
• Uniform Cartesian grids; no unstructured meshes or nesting.
• The bed changes only through setBed; sediment transport in the demos is a viewer-side illustration, not part of the verified solver.