surf

GPU kernels for breaking ocean waves, loadable through kernels: a nonlinear shallow-water solver for the surf zone and a weakly-compressible SPH solver for the overturn. The reference standards are the lake-at-rest equilibrium, held to a part in a million, and an independent torch implementation of the shallow-water right-hand side.

A swell passes through two regimes on its way to shore, and no single formulation covers both. A height field can shoal a wave, steepen it, and break it into whitewater, but eta(x, y) is a function, so it can never overturn: a barrel is a multivalued surface with air inside it. Particles can. This kernel ships the phase-resolved surf-zone hydrodynamics and the Lagrangian solver together, so a wave can be carried from the shelf to the tube, at a fifth of a millisecond per shallow-water step and 13.5 million SPH particles.

A 9-second swell driven across a 256 x 448 bed that shoals from 9 m to the beach over a submerged reef: crests steepen and refract, the reef throws a diffraction pattern downwave, and energy decays through the surf zone. On the same bed a still sea holds to 3.5e-6 m/s and 4.8e-7 m, which is what well-balancing means; the step costs 0.15 ms.

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)

version selects the release branch; trust_remote_code is required by kernels for publishers without the trusted-publisher mark.

API

Method

nsw_rhs / nsw_step solve the Saint-Venant equations with well-balanced hydrostatic reconstruction, MUSCL (minmod) reconstruction, HLL fluxes, and wet/dry runup, which shoals a swell over real bathymetry and breaks it as bores. sph_forces / build_cells run 3D weakly-compressible SPH over a uniform spatial hash: continuity with delta-SPH density diffusion, a Tait equation of state, artificial viscosity, and XSPH. Well-balancing is the property that matters most: a scheme without it spins up currents from the bed slope alone, and a still sea over a reef never settles.

Correctness

The lake-at-rest figure is measured against the published build, so it reproduces from a bare get_kernel.

Measured

Requirements 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 and a wave spills instead of tubing.
• nsw_step is explicit: dt <= 0.3 * dx / max(|u| + sqrt(g*h)).

References

Audusse et al., well-balanced hydrostatic reconstruction (2004); Harten, Lax, van Leer (HLL) fluxes; Monaghan, weakly-compressible SPH; Antuono et al., delta-SPH density diffusion.

License

Apache-2.0.