Source-linked AI summary
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
TL;DR
Large-scale water simulation lacks an efficient 2D model that is simultaneously nonlinear and dispersive, causing mismatches when coupled to nonlinear 3D solvers. The paper introduces a Hamiltonian Zakharov-based 2D solver with canonically consistent two-way coupling, reducing wave-height error and seam artifacts while substantially improving efficiency over BEM and full-domain GPU NB-FLIP.
Problem
Graphics lacks an efficient 2D wave model that simultaneously preserves nonlinear wave interactions and wavelength-dependent dispersion for coupling with nonlinear 3D flows.
Method
The method evolves the canonical Zakharov pair (η,ψ) with a Craig–Sulem DNO expansion and uses the same machinery for two-way coupling to a localized 3D Navier–Stokes solver.
Results
1.7×–2.4× lower mean wave-height error than Airy and 3–5× lower error than SWE and BEM were achieved, while the coupled pipeline showed minimal seam artifacts and ran more than 4× faster than GPU NB-FLIP.
Takeaways & Limitations
Hamiltonian two-way coupling provides a single framework for nonlinear, dispersive surface waves and localized 3D flow with improved coupling fidelity and ocean-scale efficiency.
Takeaways & Limitations
The fully nonlinear setting is not unconditionally stable, so strong 3D forcing can require reducing ε and tuning it per scene.
Abstract
from arXiv · showhide
Simulating large-scale free-surface water by coupling a localized 3D fluid solver to a cheaper 2D surface model has long faced a mismatch in wave dynamics: efficient 2D wave models used in graphics are typically either linear or non-dispersive. These models are fast, simple, and accurate for calm, small-amplitude seas, but coupling them with strongly nonlinear 3D solvers produces visible reflections and artifacts at the 2D--3D interface. We address this problem by introducing a nonlinear and dispersive 2D wave model based on the canonical Zakharov formulation. Its Hamiltonian structure, in which the surface elevation and surface potential form a canonical pair ($η$, $ψ$) governed by the wave energy, enables a canonically consistent two-way coupling scheme, allowing information to pass smoothly across the 2D--3D interface. Our 2D solver reduces mean wave-height error by 1.7--5$\times$ over SWE, BEM, and Airy baselines while running more than $10^3\times$ faster than BEM; it achieves greater nonlinear accuracy and coupling fidelity than SWE and Airy, with minor losses in speed and stability. Coupling it with a 3D Navier--Stokes solver yields a full system that suppresses visible seam artifacts across a range of experiments, including dispersion-matching and Kelvin-wake tests, and runs over 4$\times$ faster than a pure GPU NB-FLIP simulation on the same domain.
1 Introduction
The paper targets the mismatch between inexpensive 2D wave models and nonlinear 3D flow simulation by introducing a nonlinear, dispersive Eulerian model with Hamiltonian two-way coupling. This design improves wave accuracy and interface fidelity while retaining ocean-scale efficiency.
- Motivation: Large-scale free-surface simulation is computationally prohibitive with a fully three-dimensional Navier–Stokes solver, motivating localized 3D regions surrounded by cheaper 2D models.The 2D–3D decomposition is intended to propagate far-field waves efficiently while retaining 3D effects locally.
- Motivation: Existing graphics wave models typically sacrifice either nonlinearity or dispersion, limiting their ability to represent realistic open-water dynamics.Airy retains linear dispersion but discards nonlinear interactions, while SWE is non-dispersive and suited to shallow transport-dominated flows.
- Motivation: Linear 2D models coupled to nonlinear 3D solvers disagree at the interface, producing spurious reflections and visible surface artifacts.The mismatch arises because the 2D model cannot represent nonlinear effects such as crest steepening and amplitude-dependent dispersion.
- Method: The proposed Eulerian solver evolves the canonical Zakharov pair (η,ψ) under the water-wave Hamiltonian using a Craig–Sulem DNO expansion.Nonlinear corrections are evaluated with FFTs and pointwise products on a regular grid at O(N^2 log N) cost, avoiding per-step surface tracking and remeshing.
- Results: 1.7×–2.4× lower mean wave-height error than Airy and 3–5× lower error than SWE and BEM were achieved, while running more than 10^3× faster than BEM.The coupled solver produced minimal seam artifacts in dispersion-matching and Kelvin-wake tests and ran more than 4× faster than GPU NB-FLIP on the same domain.
- Method: The DNO machinery supports canonically consistent two-way coupling by lifting wave states to 3D boundary velocities and reconstituting the wave state from 3D surface information.A perturbation parameter ε continuously interpolates between linear Airy theory at ε = 0 and the nonlinear truncated Zakharov system at ε = 1.
2 Related Work
Prior work spans shallow-water, linear spectral, particle, mesh-based, and hybrid approaches, but existing coupled systems do not jointly provide nonlinear and dispersive dynamics with two-way coupling and ocean-scale efficiency.
- Surface-wave models: SWE provides efficient shallow-water simulation but lacks accurate wave dispersion for open-water applications.Later coupled methods use SWE with other solvers, but the related-work comparison identifies dispersion as a limitation.
- Surface-wave models: Airy wave theory reproduces dispersion relations but evolves linear wave dynamics, distinguishing geometric sharpening from the dynamical nonlinearity sought here.Its linearization omits nonlinear wave interactions even when surface geometry is post-processed.
- Surface-wave models: High-order spectral methods evolve nonlinear surface waves on a flat reference grid through a Taylor expansion of the Dirichlet–Neumann operator.The paper adopts this HOS expansion as the engine of its 2D wave solver.
- Alternative representations: Wave particles, packets, wavefront tracking, and surface wavelets explore alternative representations for dispersion, boundary interactions, and visual detail.These approaches are presented as related alternatives rather than as a unified nonlinear, dispersive two-way coupling framework.
- Coupled 2D–3D simulation: Prior 3D coupling methods improve efficiency or reduce reflections, but representative systems remain one-way, linear, or limited in dispersion accuracy.The FAB method is one-way, PML approaches are one-way, and recent dispersion-matching coupling still uses linear Airy theory.
- Coupled 2D–3D simulation: The proposed framework is presented as the first to combine nonlinear and dispersive surface-wave dynamics, two-way coupling, and ocean-scale efficiency in one system.Its comparison criteria include GPU-native operation, mesh-operation-free simulation, nonlinear dynamics, dispersion, and Hamiltonian evolution.
3 Hamiltonian Spectral Wave Solver
The solver evolves the canonical Zakharov surface pair with a High-Order Spectral expansion, combining nonlinear dispersive waves with FFT-based computation and coupling support. Exact linear integration with Adams–Bashforth treatment of nonlinear terms improves stability and time-step practicality, while boundary extensions and nonlinear modulation support coupled settings.
- 3.1 Canonical Zakharov formulation: The wave model treats surface elevation η and surface potential ψ as a canonical pair governed by the fluid Hamiltonian, whose energy is kinetic plus potential.The Dirichlet–Neumann operator closes the surface formulation by mapping surface potential data to the normal derivative of its harmonic extension.
- 3.2 Spectral expansion of the Dirichlet–Neumann operator: The Dirichlet–Neumann operator is approximated by a convergent Taylor expansion in η, avoiding a Laplace solve at every step and enabling FFT-based evaluation.The expansion is the High-Order Spectral method, built from Fourier multipliers and pointwise products on a regular grid.
- 3.2 Spectral expansion of the Dirichlet–Neumann operator: ε=1 gives the fully nonlinear model, ε=0 gives the linear Airy limit, and HOS-3 retains G0, εG1, and ε^2G2.The ε modulation scales nonlinear terms to improve stability in coupled simulations with steep 3D-driven surface features.
- Boundary treatment: Even extensions on doubled grids support reflecting and replicating boundaries, enforcing zero discrete normal derivatives for reflecting walls while increasing FFT cost only by a constant factor.After each step, the solution is cropped back to the original grid.
- 3.4 Time integration: exact-linear Adams–Bashforth: The exact-linear Adams–Bashforth scheme integrates linear dispersion exactly per Fourier mode while advancing only the nonlinear forcing explicitly.This removes high-wavenumber linear stiffness and reduces the time-step restriction to a milder CFL condition on the nonlinear terms.
- 3.4 Time integration: exact-linear Adams–Bashforth: The same exponential-integrator strategy gives closed-form treatment to the position-dependent relaxation source used for coupling the wave surface toward a target state.The source is linear in η−ηtarget and is composed with the wave step.
4 2D-to-3D Coupling: Bulk Velocity from the Canonical Surface State
The 2D-to-3D coupling reconstructs bulk velocities from the canonical surface state using Craig–Sulem expansion and finite-depth harmonic kernels, then supplies them to the 3D boundary layer.
- Bulk velocity coupling: The localized 3D solver receives boundary velocities reconstructed from the canonical surface pair (η, ψ) of the surrounding 2D wave model.The 3D region is embedded within the open-ocean wave domain, with its lateral and bottom boundaries supplied by the wave state.
- FAB boundary layer: FAB cells exclude the 3D pressure projection and receive Dirichlet data from the wave model, including the level set ϕ(x,y)=y−η(x,t).The FAB band is four cells wide in all runs.
- Potential reconstruction: Craig–Sulem reconstruction extends the surface potential ψ to the flat reference plane before harmonic continuation through the finite-depth fluid strip.The expansion is truncated at first order, Φ|y=0≈ψ−εηG0ψ, to remain robust under arbitrary 3D forcing.
- Finite-depth kernels: Finite-depth kernels satisfy the surface and bottom conditions, recovering G0 at the surface and zero vertical velocity at the rigid bottom.In deep water, both kernels reduce to exponential attenuation, generalizing linear deep-water coupling to finite depth.
- Velocity transfer: The reconstructed velocity field is sampled at multiple depth layers, transformed by inverse FFT, and trilinearly interpolated onto the 3D FAB boundary cells.The cost is Ld FFT pairs per substep on the shared surface grid.
5 3D-to-2D Coupling: Implicit Relaxation in Canonical Variables
The 3D-to-2D direction transfers the 3D surface elevation through a spatially graded relaxation source, while reconstructing the canonical potential to preserve consistency.
- 3D-to-2D transfer: The 3D solver’s surface elevation is transferred back to the 2D wave field through a relaxation zone rather than by overwriting the wave state.This allows hull depressions, splashes, and other 3D features to enter the surrounding wave field smoothly.
- Target construction: The target elevation is extracted by a column-wise level-set scan, then padded and Laplace-smoothed to suppress staircase noise.The surface cell is identified after three consecutive air/solid cells are detected.
- Spatial relaxation: The relaxation rate is zero outside the box, rises smoothly across the transition band, and reaches σmax inside the embedding region.Because the rate depends on position rather than the wave field, the source is linear and exactly integrable.
- Canonical consistency: Only η is relaxed; the wave’s surface velocity is preserved and ψ is recomputed through an inverse DNO solve so the canonical state remains self-consistent.An independent velocity target risks introducing spurious reflections and over-constraining the canonical pair.
- Stable integration: The exact exponential relaxation is unconditionally stable and monotonically convergent, avoiding the explicit restriction σmaxΔt<2.Operator splitting applies the Hamiltonian wave step first and the pointwise relaxation second.
6 Implementation Details
The implementation uses GPU spectral operations and post-processing while acknowledging that demanding coupled scenes can retain visible boundary reflections.
- GPU implementation: The solver uses Warp for GPU acceleration and PyTorch FFT and padding operations for the 2D spectral computation.It supports standard ocean-spectrum initialization, including TMA spectra with Donelan–Banner directional spreading.
- Rendering: Whitewater and fine-scale surface detail are added as post-processes rather than included in the quantitative simulation comparisons.Whitewater uses Houdini’s solver, while render-time ripples are disabled for quantitative evaluation.
- Residual artifacts: In the heaviest-seas battleship scene, wavefronts still reflect from the rectangular 3D-box boundary and propagate into the 2D domain.Render-time spectrum detail and foam can obscure these residual artifacts in production settings.
- Spectral robustness: The 3/2-rule de-aliases nonlinear products by padding spectral factors, multiplying in physical space, and truncating back to the resolved grid.A smooth taper additionally limits high-wavenumber noise amplification from nonlinear corrections and coupling boundaries.
- Velocity filtering: Spectral tapering suppresses unreliable high-wavenumber forcing before the reconstructed velocity is passed to the 3D pressure solve.The taper preserves the lower half of the resolved spectrum while decaying through the upper half.
7 Time Integration
The coupled algorithm advances 3D motion and the canonical 2D wave state in a coordinated substep loop, with optional 2D sub-cycling and bidirectional state transfer.
- Time stepping: The 3D solver sets the outer timestep from advection and solid-motion CFL constraints, while the 2D solver may use the same step or sub-cycle with a smaller one.Most scenes use Δt2D=Δt3D.
- 3D update: Each cycle advances 3D transport and solid motion, then extracts a smoothed target elevation from the tentative 3D level set.The target scan uses three consecutive air/solid cells before reading the surface.
- 2D update: The 2D canonical state is advanced under the Zakharov Hamiltonian using an integrating-factor AB2 scheme before applying relaxation toward the 3D target.The post-advance state is mapped forward through the DNO before relaxation.
- State consistency: After relaxation, the inverse DNO map reconstructs ψ from the updated η and preserved vertical velocity, maintaining a consistent canonical pair.The resulting height is embedded into the 3D FAB ring through ϕ=y−η.
- Coupling loop: The algorithm reconstructs depth-layered bulk velocities from the 2D state, interpolates them into the FAB ring, and advances the 3D solver with those boundary conditions.This closes the two-way coupling loop between the surface model and localized NB-FLIP region.
8 Validation
Validation shows that the nonlinear, dispersive wave model improves standalone accuracy over SWE, BEM, and Airy baselines, preserves nonlinear wave behavior, and couples with 3D flows with fewer interface artifacts. The method remains stable across demanding scenes but requires parameter tuning and can exhibit reflections near-breaking forcing.
- Wave propagation and nonlinear accuracy: 1.7×–5× lower mean wave-height error than Airy, SWE, and BEM is achieved while remaining more than 10^3× faster than BEM.The strongest linear baseline is improved by 1.7× in deep water and 2.4× in shallow water; SWE and BEM improve by roughly 3–5×.
- Wave propagation and nonlinear accuracy: Order 2 and Order 3 HOS reduce wave-height error monotonically as ε increases toward the fully nonlinear model.The decrease is observed against a 3D Eulerian reference in shallow-water propagation and against analytical Stokes solutions across tested initial conditions.
- Wave propagation and nonlinear accuracy: Fully nonlinear HOS-3 remains phase-locked to the analytical Stokes wave, whereas the linear solver progressively drifts out of phase.For steeper waves, Order 3 outperforms Order 2; for milder waves, the two truncations are nearly indistinguishable.
- 2D–3D coupling validation: The coupled model maintains symmetric wave propagation without artificial reflections at the 2D–3D interface, unlike models that fail to capture dispersive dynamics.The comparison includes a wave packet crossing the interface and a 2D region where SWE produces incorrect propagation.
- Stress tests and ablations: Disabling the low-pass filter quickly destabilizes Crown Splash, while first-order integration increases shallow-water mean/max error from 0.011501/0.051843 to 0.012752/0.053980.The ablation isolates stability and accuracy contributions from the filter and IF-AB2 integrator.
- 2D–3D coupling validation: The full method reproduces feather-like Kelvin wakes with minimal interface reflection and runs more than 4× faster than pure GPU NB-FLIP.The reported per-substep times are 85.0 ms for the coupled method and 363 ms for NB-FLIP.
- Stress tests and ablations: Near-breaking heavy-wave scenes require ε=0.2 and can still produce reflected wavefronts from the rectangular 3D-box boundary.The restricted nonlinearity cannot fully absorb strongly nonlinear content generated by the 3D solver.
9 Conclusions and Limitations
The paper concludes that a Hamiltonian nonlinear and dispersive 2D model improves coupling fidelity in strongly nonlinear scenes at the cost of speed, stability, and representational scope. The authors identify scene-dependent nonlinearity tuning, truncation overhead, and inability to represent breaking waves as key limitations.
- Conclusions: The solver trades modest speed and stability losses for greater nonlinear accuracy and coupling fidelity than prior linear or non-dispersive 2D models.The method is intended for strongly nonlinear scenes, while simpler models remain preferable for gentle, nearly linear surfaces.
- Conclusions: The canonical Zakharov pair is advanced under the exact water-wave Hamiltonian using a Craig–Sulem expansion with O(N^2 log N) cost.The nonlinear corrections reduce to FFTs and pointwise products on a regular grid.
- Limitations: Fully nonlinear ε=1 is not unconditionally stable, so strong 3D forcing may require reducing ε to maintain stability.Experiments use ε=0.3 for the submarine and ε=0.2 for the battleship and pond.
- Limitations: Order 3 is roughly 2× slower than Order 2, while Order 2 is roughly 5× more expensive than linear-wave solvers.Order 3 provides only marginal accuracy gains in the steepest regimes, and users must tune ε per scene.
- Limitations: The 2D height-field model cannot represent breaking waves once the surface overturns.At overturning, the single-valued surface elevation representation is no longer valid.