Simulating Shallow Water Wave Equations in GLSL

Simulating shallow water wave propagation on 2D heightfield grids in GLSL is most effectively achieved using the Explicit Finite Difference Method (FDM) applied to the discrete 2D wave equation or the full Saint-Venant shallow water equations. By utilizing GPU framebuffer ping-ponging, fragment shaders compute height and velocity changes across a grid in real time, making it the standard approach for interactive water surfaces in games and real-time graphics.

The Algorithm: Explicit Finite Difference Method (FDM)

The algorithm discretizes space and time on a uniform 2D grid. While the full Saint-Venant equations track both fluid height (\(h\)) and horizontal velocities (\(u, v\)), real-time heightfield graphics commonly use the linearized 2D wave equation approximation:

\[\frac{\partial^2 h}{\partial t^2} = c^2 \left( \frac{\partial^2 h}{\partial x^2} + \frac{\partial^2 h}{\partial y^2} \right)\]

Where:

  • \(h\) is the fluid surface height.
  • \(c\) is the wave propagation speed, derived from gravity and depth (\(c = \sqrt{gh}\)).

Numerical Discretization

To solve this in a GLSL fragment shader, continuous derivatives are replaced with finite central differences. Using a 5-point discrete Laplacian stencil over adjacent grid cells (texels), the height at the next time step (\(t + \Delta t\)) is derived as:

\[h^{t+\Delta t}_{x, y} = 2h^t_{x, y} - h^{t-\Delta t}_{x, y} + \left(\frac{c \Delta t}{\Delta x}\right)^2 \left( h^t_{x+1, y} + h^t_{x-1, y} + h^t_{x, y+1} + h^t_{x, y-1} - 4h^t_{x, y} \right)\]

To ensure physical realism and numerical stability, a damping factor (\(d \in [0.98, 1.0)\)) is applied to simulate viscous dissipation.

GLSL Ping-Pong Implementation

Because GPU shaders cannot write to and read from the same texture simultaneously, the simulation employs a ping-pong rendering technique using two (or three) floating-point Framebuffer Objects (FBOs):

  1. State Storage: The current state texture stores the current height (\(h^t\)) and previous height (\(h^{t-\Delta t}\)) in separate color channels (e.g., Red and Green).
  2. Fragment Shader Execution: A full-screen quad is rendered. For each pixel:
    • Sample the current height of the north, south, east, and west neighbors.
    • Calculate the discrete Laplacian.
    • Compute the new height (\(h^{t+\Delta t}\)).
    • Apply damping: \(h^{t+\Delta t} = h^{t+\Delta t} \times d\).
    • Output the new height to the target buffer's Red channel, and write the current height (\(h^t\)) to the Green channel.
  3. Swap Buffers: Swap the input and output textures for the next frame.

Stability Requirement: CFL Condition

For the explicit finite difference scheme to remain stable without exploding into numerical artifacts, it must satisfy the Courant-Friedrichs-Lewy (CFL) condition:

\[C = \frac{c \Delta t}{\Delta x} \le \frac{1}{\sqrt{2}} \approx 0.707\]

If the wave speed \(c\) or the frame time step \(\Delta t\) is too large relative to the grid spacing \(\Delta x\), the simulation will destabilize. In GLSL implementations, clamping \(\Delta t\) or running multiple simulation sub-steps per frame ensures ongoing stability.