In Wavebreakers, the water is the game. An arcade jetski racer depends on exciting, fun, and accurate water interactions. Because of this, the water gets the lion’s share of the performance budget. I’ve allocated 12ms out of a 16.7ms frame time, assuming 60fps target — everything else is secondary and must fit into the remaining 4ms.
A key architectural decision drives the simulation design: water must be CPU-analytic and deterministic. And, as much as possible, nothing gameplay relevant should depend on prior state.
Every buoyancy probe resolves against a deterministic maths function on the CPU, dependent only on (x, z) world coordinates, and t (time). The visuals follow the physics; layers of visual effects are added on top (to be covered in a future post), but the simulation is always the source of truth.
The failure of the vertical sine wave
Intuitively, most developers (myself included) start by imagining an ocean as a field of sine waves. Sample a height as , displace each mesh vertex vertically by that amount, and you’ve got an ocean. Normals are trivial to calculate, so shader math becomes cheap. Indeed, my first attempt at an ocean game took this approach, but immediately I realised it looks rather underwhelming. A vertical sine wave rolls like a rubber sheet. Every crest is a rounded hump and every trough its exact mirror, because a vertical-only sine displacement is symmetric. Humans are particularly tallented at spotting repeating patterns, and it immediately reads as repetitive, fake, and unnatural.
Real ocean swell is not symmetric — crests pinch into sharp ridges, troughs stretch and flatten — no amount of amplitude tuning or frequency layering creates that shape from terms that only move in one dimension.
If you’ve ever been unlucky enough to study signals engieering, you’re likly thinking a Fourier series can approximate any periodic function, so the answer is just to add more waves. There are two issues with this approach.
First, the performance cost scales quickly. This approach naturally lends itself to running in a GPU vertex shader, where even high-resolution meshes and large Fourier series can be computed in parallel extremly cheaply on modern hardware. The problem with that approach is the true surface now exists only on the GPU, where the CPU based physics can’t access it. As far as the CPU is concerned, the ocean surface is still a flat mesh. GPU readbacks are possible, but they’re extremly expensive and, crucially, non-deterministic due to the parallel nature of GPUs. So either the simulation becomes decoupled from the visuals, or the physics needs to run on the GPU too — a path that is not viable for a multiplayer game where every client must evaluate the same surface at the same time (not to mention an order of magnitude more complex to implement and debug). The only practical way to get a physically readable surface, deterministically and in lock-step, is to compute it on the CPU - and that means the maths must be cheap enough to run on the CPU for every hull probe, every tick, for every player in an online race.
The second issue is that, no matter how many sine waves we add, the surface remains symmetric. The up-and-down-only displacement is a fundamental limitation of the sine wave approach.
Enter, the Gerstner wave (aka the trochoidal function)
One way to break that symmetry is to add a horizontal component to the displacement. A Gerstner wave displaces the surface both vertically and horizontally, which pinches crests and stretches troughs, and gives the whole surface a sense of surging and rolling. There are other advantages, which we’ll touch on later (useful for calculating where to add foam and spray effects, for example), but the key advantage is the addition of the horizontal term to break the symmetry.
Not being the best study when it came to maths, I did find it significantly more difficult to get my head around at first, so let’s break it down.
Rolling in circles
You can think of a Gerstner wave as a vertical sine wave with a horizontal cosine term added. And if you remember your trigonometry, that’s literally working backwards to a circle (sine being the vertical component of a swept unit circle, and cosine being the horizontal).
The vertical term moves the surface up and down, while the horizontal term nudges the surface side to side. The result is that the surface point (say, a mesh vertex in our case) traces a circle, rather than just moving up and down. When calculated over neighbouring points on the mesh, the horizontal term pinches the crests and stretches the troughs, while the vertical component controls wave height. This does a pretty good job of approximating real ocean swell, without using any complicated or expensive fluid simulation techniques.
By altering the ratio between the horizontal and vertical terms we can control the steepness of the wave. And, much like the approach of adding sines, we can add multiple Gerstner waves together, altering parameters such as wavelength, phase, direction, amplitude, and steepness, to create a complex sea state; remaining cheap enough to run on the CPU (or even in a web browser - the hero banner on the homepage of this blog is actually a live Gerstner simulation).
We get better looking waves than a sum of sines, and we have more individual knobs to tune to create varied and interesting sea conditions. Win-win.
Constructing the sum of components
So our sea state is now a sum of Gerstner components — each a directional wave carrying its own parameters, all added together — evaluated as a pure function of position and game time.
Because I’m trying to emulate a realistic (if stylized) ocean, wave speed isn’t treated as a free parameter, because long ocean waves travel faster than short ones. To make it look convincing I derive each component’s speed from its wavelength.
Stacking these waves is vector addition. A surface point’s total offset is the sum of each orbital vector. Or to put it another way, the end point of one wave’s orbit is the origin point of the next — and the end of the chain is the water surface itself.
For Wavebreakers, I’ve chosen to use twelve waves. There is no mathematical basis for this number, it’s simply the point beyond which I felt there were diminishing returns, and is low enough that even a modest system can keep up. It may well change as I continue to develop the simulation.
The spectrum I use is most dominant in the long wavelengths, with amplitude falling as roughly λ^1.7, and wavelengths roughly following a log spread. I add some jitter, and speeds are then speeds are tuned so the components don’t lock into step with each other.
The exact values are mostly arbitrary, but the offsets are carefully designed to avoid artifacts which emerge when multiple near-equal components tile into a regular pattern (if you’ve ever heard two musical notes played, where one is slightly out of tune with the other, this is the beating pattern you hear).
Decorrelated phases ensure patterns never develop, long-wave dominance ensures it resembles a natural ocean surface.
The vector displacement formula
For one component, every point point on the flat plane maps to a displaced surface point:
- is amplitude
- the wave’s unit heading
- the angular frequency
- a per-component phase offset (this is used to ensure every component doesn’t stack at the world origin, causing a visible seam)
- the steepness coefficient
The actual code loop is summed over the spectrum of waves. The result is a single function that takes a surface position and a game time, and returns the displaced surface position:
# sim/water/water_surface.gd — the one displacement, per component
for i in waves:
var w := _a[i] # dir_x, dir_z, k, omega
var ab := _b[i] # amplitude, Q, phase, breath
var theta := w.z * (w.x * px + w.y * pz) - w.w * t + ab.z
var c := cos(theta)
var qa := ab.y * amp # Q · A — Q already normalised (see below)
dx += qa * w.x * c # horizontal push, along the heading
dz += qa * w.y * c
dy += amp * sin(theta) # vertical riseThere are some additional details in the implementation, such as a “breath” term to slightly modulate the amplitude of each component over time, normalisation of the steepness to prevent the horizontal push from exceeding the self-intersection limit (the point at which the wave has been pinched so far it folds over itself), and attenuation of each component’s amplitude and steepness as it approaches the shore.
We’ll cover these aspects in a later article.
Where physics breaks: the self-intersection limit
As useful as the Gerstner wave is, the horizontal squeeze does bring some consequences we have to deal with. The Jacobian of the horizontal map (a way to tell if neighbouring points keep their left-to-right order) works out to:
.
- is the wavenumber
- the amplitude
- the steepness.
If steepness stays low, the result is positive everywhere, meaning the surface is single-valued — only one height for every (x, z) position — and there’s a definitive answer to the question “how high is the water here?”.
Once Q goes past this limit and the Jacobian goes negative, neighbouring points start to
cross (the surface folds back through itself), and the surface stops being a function of (x, z).
In this case the surface is now multi-valued, and the simulation can no longer answer “how high is the water here?”, because there are multiple answers. This steepness limit is a function of the sum of over the spectrum, and is a hard limit: if any component’s , the surface will fold.
For a single wave that limit is exactly , which is why the
steepness dial is normalised by and authored as a fraction of the fold
(0.6 = 60% of the way there), rather than an absolute number.
When we start adding multiple components there’s a chance that every crest might align at a single point. I use a weighted safety factor to clamp the sum slightly just slightly past the fold limit.
The surface is only mostly prevented from folding, because we can actually use the fold to our advantage, as long as we don’t let it visually break the illusion of a continuous surface.
By restricting the fold to a small fraction of the surface at only the steepest crests, we can use the negative Jacobian to analytically determine where foam and spray should appear (eg, where the tips of the waves are breaking), and use visual effects to hide the fold itself. In this way, the foam and spray are accurate representations of a real physical phenomenon (a wave breaking), and the limitation now becomes a feature.
Since a folded surface doesn’t render well, I use J < 0 to paint crest foam and generate spray particles - hiding the fold and selling the illusion of rough, churned sea.
Working backwards: finding the surface height at a given (x, z)
Adding the horizontal component does make it more comlicated to get a surface height value at a specific (x, z) coordinate, because the horizontal push means the surface is
no longer a simple function of (x, z).
In the forwards direction (calculating the final surface point position from a starting position), we just need to plug in the (x, z)
starting coordinate and the formula will give us back the displaced surface position. This is what’s used in the vertex shader to warp the flat ocean surface mesh into the correct shape.
There are times we need to be able to work backwards though, such as calculating final surface height at a specific (x, z) coordinate - most notably for buoyancy calculation in the physics sim.
With the sine wave approach, vertical displacement is simply the sum of sines, and is trivial to invert: given an (x, z) you can directly compute the height, because there is no horizontal translation, so this point is the parameter passed into the height function.
With a Gerstner wave, the horizontal displacement means that the height of the final surface that corresponds to a given (x, z) coordinate does not originate at that same (x, z), and we need a way of calculating the correct input coordinate.
To find the surface height, we have to solve for the starting point on the flat plane (which changes every time step). This is achieved by a method of successive approximation, taking multiple iterations to zero in on the correct values (Newton-Raphson method).
This extra computation can stack up, as buoyancy for the Jetskis is calculated by reading six probes per hull, and each probe (and any other items calculating buoyancy) needs to find the correct surface height at its (x, z) location. In practice this is trivial on modern hardware, and the performance is well within budget even for a multiplayer game.
Achieved determinism at scale
This approach gives us a highly tunable system, but keeps the entire surface as a pure function of (x, z, t). There is no stored state, no
per-frame integration, and nothing to potentially accumulate error by building on the previous state.
Every consumer calls the same function and gets the same answer: server simulation, client prediction, reconciliation replay, buoyancy probes, etc. This shields us from any network divergence issues.
Because the surface is evaluated rather than simulated, there is no way for it to
diverge,in the way a stateful physics sim can drift.
The host’s clock is the single source of the t everyone evaluates, and the host’s spectrum is the single source of the wave parameters everyone uses. The result is a deterministic, lock-step simulation that can be run on any number of clients, with no divergence, no need for reconciliation, and no heavy state data to send over the network.
Next up, we’ll dive into some shader code, and see how this foundation is used to generate foam, spray, and caustics, and how the layering of these visual effects creates a convincing ocean.