Hamiltonian Two-Way Coupling of Nonlinear Waves and 3D Flows
Sinan Wang, Ruicheng Wang, Taiyuan Zhang, Fan Feng, Jinjin He, Yuchen Sun, Zhiqi Li, Bo Zhu
SIGGRAPH Asia 2026)
cs.GR, math.NA, physics.flu-dyn
2026-08-26
Zakharov 2D waves two-way coupled to local 3D NS cut mean height error 1.7–2.4× vs Airy and run 4× faster than full-domain GPU NB-FLIP.
A full-domain 3D incompressible Navier-Stokes solve blows up in memory and compute once the water spans hundreds of meters. Graphics practice puts the expensive 3D solver in a small box around a hull or splash, then lets a cheap 2D surface model carry far-field waves. That split only works if both sides agree on the surface state at the seam.
The 2D models graphics actually uses give up one of two things. Airy theory keeps the right dispersion, so different wavelengths travel at different speeds, but it linearizes the free surface and drops wave-wave interaction that grows with steepness. The shallow-water equations keep nonlinear advection and collapse dispersion to ω = √(gh) k, which is fine in a tank and wrong in open water. Wave particles, packets, and front trackers ride linear rays and inherit the same linearization. A 3D box, by construction, sends nonlinear waves to its boundary: crests steepen, phase speed depends on amplitude. A linear 2D model cannot hold that state, so the seam reflects. Boundary-element methods recover nonlinear potential flow, then pay O(N³) dense interactions plus remeshing.
Georgia Tech and Dartmouth replace the far field with Zakharov's canonical pair: surface elevation η and surface velocity potential ψ. Total wave energy is the Hamiltonian, so (η, ψ) evolve like position and momentum. The one nontrivial object is the Dirichlet-Neumann operator G(η): given ψ on the moving surface, it returns the normal derivative of the harmonic extension of the potential into the fluid. With G in hand, the 3D Laplace problem collapses to two fields on a flat horizontal grid.
Evaluating G directly would solve Laplace in an irregular body of water at every step. They use the Craig-Sulem High-Order Spectral expansion instead, writing G as a series of homogeneous powers of η. Each term is Fourier multipliers plus pointwise products, so the cost is O(N² log N) FFTs. Production runs truncate at order 2 (HOS-2) or order 3 (HOS-3). A scalar ε ∈ [0, 1] then scales every nonlinear correction: ε = 0 is linear Airy, ε = 1 is the truncated nonlinear Zakharov system. When 3D forcing pushes the surface toward breaking, they turn ε down for stability.
Time integration rotates the linear part exactly in Fourier space at each mode's dispersive frequency, then advances the nonlinear residual with Adams-Bashforth-2. Linear waves are therefore non-dissipative, and the time step is no longer pinned by the stiffest wavenumber. Reflecting walls use even extension on a doubled grid; nonlinear products use 3/2-rule de-aliasing.
Both directions of coupling reuse the same DNO:
Standalone 2D test: a 20 m × 20 m tank, 1 m shallow and 20 m deep, a Gaussian bump evolved for 2000 frames against a high-resolution 3D Eulerian level-set reference. Time-averaged absolute height error (shallow / deep):
| Method | Shallow mean | Deep mean |
| SWE | 0.0482 | 0.0962 |
| BEM | 0.0517 | 0.1705 |
| Airy DK | 0.0277 | 0.0545 |
| HOS-2 | 0.0136 | 0.0328 |
| HOS-3 | 0.0115 | 0.0342 |
Against the strongest linear baseline, Airy DK, mean error drops about 2.4× in shallow water (HOS-3) and 1.7× in deep water (HOS-2). Against SWE and BEM the factor is about 4× shallow and 3 to 5× deep. On an RTX 4090, HOS-2 takes 6.2–6.4 ms per step and HOS-3 about 13.5 ms. Airy DK is about 1.1 ms; BEM is 8.7 s deep and 26.0 s shallow. That is three orders of magnitude faster than BEM, and roughly 5× slower than the linear solvers.
On Stokes wavetrains, the ε = 0 linear solver lacks the amplitude-dependent frequency shift and drifts out of phase. HOS-3 at ε = 1 stays locked to the analytical second-order Stokes profile. In the dispersion-matching test, the SWE coupling of Chentanez et al. (2015) is left-right asymmetric. Schreck and Wojtan (2022), even after retuning their dispersion factor to d = 0.8, is still less consistent across the seam. FAB + Airy DK is symmetric but produces deeper waves inside the 3D box. A single-boat Kelvin wake is the wrong shape under SWE, and several Airy or hybrid couplings leave spurious stern reflections. This pipeline carries the wake through the seam at 85.0 ms per substep, versus 363 ms for a full-domain GPU NB-FLIP of the same visible water, more than 4× faster. Larger shots (seaplane, battleship, surfacing submarine) keep the seam mostly clean. ε is 0.3 on the submarine, 0.2 on the battleship and pond, and 1.0 on most other scenes.
Prior 2D-3D couplings drop nonlinearity, drop dispersion, or pay BEM's mesh cost. This pipeline stacks all four with structured-grid FFTs. For hull wakes and splash radiation, the 3D box handles breaking and solids while the 2D grid carries dispersive, nonlinear waves without an NB-FLIP bill for the whole sea. ε is a continuous knob from linear Airy to truncated Zakharov, easier to tune than swapping solvers.
Calm, nearly linear far fields still want Airy or a Tessendorf spectrum: cheaper and more stable. HOS-2 is enough for almost every scene they ran; HOS-3 adds only a little accuracy in the steepest regimes.
ε = 1 is not unconditionally stable. Strong 3D forcing near the breaking threshold lets higher-order DNO terms blow up, so ε has to be turned down per scene. Each extra truncation order roughly doubles the wave-step cost. The 2D model is a single-valued height field, so overturning breakers are out of representation. The whole 2D grid advances every substep, so the domain cannot grow indefinitely the way a Tessendorf spectrum can. Residual seam artifacts remain on the battleship in heavy seas; they hide them at render time with extra spectral displacement and foam.
Ablations confirm the coupling is intentionally asymmetric: pinning η is stable, pinning ψ is not, and disabling the spectral filter kills Crown Splash quickly. The paper never reports a numeric seam-error metric; dispersion matching and the boat comparisons are visual. The BEM timings at 8.7–26 s also depend on implementation and meshing. ε is still a manual knob, with no adaptive policy.