surf

GPU kernels for breaking ocean waves, loadable through kernels. Two solvers, covering the two regimes a swell passes through on its way to shore.

nsw_rhs / nsw_step — nonlinear shallow-water (Saint-Venant). Well-balanced hydrostatic reconstruction, MUSCL (minmod) reconstruction, HLL fluxes, wet/dry runup. Shoals a swell over real bathymetry and breaks it as bores: the surf zone, the whitewater, the runup on the sand.

sph_forces / build_cells — 3D weakly-compressible SPH over a uniform spatial hash. Continuity with delta-SPH density diffusion, Tait equation of state, artificial viscosity, XSPH. Lagrangian, so the free surface is free to overturn.

Why both

A shallow-water model (like a spectral one) carries a single-valued surface. eta(x, y) is a function, so it can shoal a wave, steepen it, and break it into whitewater, but it can never overturn: a barrel is a multivalued surface with air inside it, and no height field can represent that. That regime needs particles. surf ships the phase-resolved surf-zone hydrodynamics and the Lagrangian solver for the overturn, so a wave can be carried from the shelf to the tube.

Validation

The lake-at-rest and free-fall figures are measured against this published build (v1, sm_86), not a local reference, so they reproduce from a bare get_kernel.

Well-balancing is the one that matters most: a scheme that is not well-balanced spins up currents out of the bed slope alone, and a still sea over a reef never settles. This one holds still to a part in a million.

Usage

import torch
from kernels import get_kernel

surf = get_kernel("phanerozoic/surf", version=1, trust_remote_code=True)

# surf zone: shoal and break a swell over bathymetry (b = bed elevation)
for _ in range(nsteps):
    h, hu, hv = surf.nsw_step(h, hu, hv, b, dx, dy, dt, cf=0.004)
eta = h + b                      # free surface; h = 0 is dry sand

# overturning: WCSPH forces (pos/vel [N,3], ptype 0=fluid 1=wall)
perm, cs, ce = surf.build_cells(pos, origin, csz, grid)
pos, vel, rho, ptype = pos[perm], vel[perm], rho[perm], ptype[perm]
drho, acc, xsph = surf.sph_forces(pos, vel, rho, ptype, cs, ce,
                                  grid, origin, csz, h=hsm, mass=mass, c0=c0)

Scope and limits

• float32, CUDA, compute capability 8.0+, Linux x86_64 (the kernel builder emits no Windows variants).
• Shallow-water fields are [Ny, Nx]; x is cross-shore (zero-gradient), y is longshore (periodic). b is bed elevation, eta = h + b, h is depth.
• NSW is depth-averaged and non-dispersive: right for the surf zone, wrong for deep-water swell propagation over long distances.
• SPH is weakly compressible; c0 should be ~10x the fastest expected flow. Resolving a plunging jet needs roughly H/dp >= 25; coarser than that and a wave spills instead of tubing.
• nsw_step is explicit: dt <= 0.3 * dx / max(|u| + sqrt(g*h)).

License

Apache-2.0.