Ocean Simulation

This project focus on simulating the ocean in real time using the Inverse Fast Fourier Transform technique. This is implemented in Unity URP, utilizing the GPU compute shader and Stockham FFT for GPU for efficiency. The project is largely based on Jerry Tessendorf's popular SIGGRAPH Course Simulating Ocean Water.^[1]

Related Works

Jerry Tessendorf^[1] and Kass and Miller^[2] were among the first to simulate water surfaces by solving the wave equation with internal boundaries on a 2D height field.

Horvath ^[3]described the practical application of several empirically-based, directional ocean wave spectra, extending the spectrum functions in Tessendorf's paper.

Chentanez and Müller ^[4] further specialized shallow water solver that can handle arbitrary underlying terrain slopes, arbitrary water depth and supports wet-dry regions tracking.

Statistical Wave Models with Inverse Fourier Transform

Figure 1: Nice diagram for fourier series from Duke Institute for Brain Sciences Methods Meetings.

Build a complex wave using simple sine waves

Statistical models from the literature of oceanography tend to represent the wave height field as a sum of sine and cosine waves. This is analogous to the concept of 1D inverse Fourier transform.

Figure 1: Representation of wave height field as a summation of sine and cosine waves

By applying this representation, real-life measurement of spectra is connected with the model, making the visualization of the ocean more realistic. The FFT-based representation of a wave height field expresses the wave height $h(\mathbf{x}, t)$ at the horizontal position $\mathbf{x} = (x, z)$ as the sum of sinusoids with complex, time-dependent amplitudes:

h(\mathbf{x}, t) = \sum_{\mathbf{k}} \tilde{h}(\mathbf{k}, t) \exp(i \mathbf{k} \cdot \mathbf{x})

where t is the time and

\mathbf{k}

is a two-dimensional wave vector with components

\mathbf{k} = (k_x, k_z)

, describing the wave that propogating towards unit vector

\mathbf{k}

with wavenumber

|\mathbf{k}|

.
The wave vector's component along x and z must be integer multiple due to the orthogonality of the fourier basis.

k_x = 2πn/ L_x, k_z = 2πm /L_z

, where n and m are integers with bounds

-N/2 \leq n \lt N /2, -M/2 \leq m \lt M/2

.
The FFT process generates the height field at discrete points

\mathbf{x} = (nL_x / N, mL_z / M)

Phillips Spectrum - How much portion of each wave we add into our ocean

After reaching this expression, we can combine the energy at each wave vector $\mathbf{k}$ with measurements from real life. There are many spectra, measured using a number of wave-buoy, photographic, or radar measurements of the ocean surface. One of the simplest and earliest one is called Phillips Spectrum. It is measuring the intensity at each wavenumbers $\mathbf{k}$ , narrating the fact that the wave height amplitudes $\tilde{h}(\mathbf{k}, t)$ in the real world are nearly statistically stationary, independent, Gaussian fluctuations with a spatial spectrum denoted by:

P(\mathbf{k}) = A \frac{e^{-1/(k L)^2}}{k^4} |\hat{\mathbf{k}} \cdot \hat{\mathbf{w}}|^2

L = \frac{V^2}{g}

The meaning of those parameters are as follows:

$\mathbf{k}$ - Wave vector.
$\hat{\mathbf{k}}$ - Wave direction, normalized wave vector.
$k$ - Wave number, norm of wave vector. Same quantity as mentioned before.
$A$ - A numerical constant.
$\hat{\mathbf{w}}$ - direction of the wind.
$L$ - Largest possible waves length arising from a wind of speed V.
$g$ - Gravitational constant.
$V$ - Wind speed.
$L$ - Largest possible waves length arising from a wind of speed V.

The intuition behind the equation is quite clear and remember it is a empirical model:

ExplainationOfPhillipsSpectrum — Figure 2. Explaination of Phillips Spectrum.

Although its simplicity, this model has poor convergence properties at high values of the wavenumber |k|. A simple fix is to modify the spectrum by the multiplicative factor: $\exp(-k^2L^2)$ .

Figure 3. Phillips spectrum with wind speed L = 3 and directional spreading term D = 0.06. One observation from the function with different parameters: the spectra suppresses both the low and high wavenumber waves. (Which is a quite common phenomenon in different aspects in real life.)

Get It Moving! -Dispersion Relation

Since we have a way to represent complex wave and a way to get the proportion of each wave. The next motivation is to make these wave move. Following is how we achieve it.
After we get the spectra for waves, we need to generate the initial waves based on the spectra. Based on the assumption that the height amplitudes in the real world are nearly stationary and independent Gaussian fluctuations, we can generate the power at each wave vector under the assumption of Phillips Spectrum using the equation:

\tilde{h}_0(\mathbf{k}) = \frac{1}{\sqrt{2}} \left( \xi_r + i\xi_i \right) \sqrt{P_h(\mathbf{k})}

where $\xi_r$ and $\xi_i$ are ordinary independent draws from a Gaussian random number generator, with mean 0 and standard deviation 1. This equation keeps the expectation of the amplitudes equal to $P_h(\mathbf{k})$ while introducing the properties of Gaussian from the real-world observation.

The only thing left is how to animate the waves.
There are 2 key concepts to make it possible. First, from the dispersion relation derived from the Navier-Stokes equations, we know the relation between the wavenumber $k$ and the wave temporal frequency $\omega(k)$ :

\omega^2 = gk

This means larger wave numbers leads to higher temporal frequency (angular frequency - how frequently a point on the water surface moves up and down over time). However, the derivation of the relation is quite complicated involving linearization of the differential equations, which I don't understand :). (Mark date: 11/4/2025).

The second key concept is the temporal phase evolution of each wave component. This is done by applying phase shift to the wavevector $\mathbf{k}$ by its temporal frequency, and applying the reverse phase shift to the wavevector in the opposite direction. This also ensure a real-valued surface height. (Phase shift is rotation in complex domain, which is done by multiply a complex number with exponential complex number)

\tilde{h}(\mathbf{k}, t) = \tilde{h}_0(\mathbf{k}) \exp\{i \omega(k) t\} + \tilde{h}_0^*(- \mathbf{k}) \exp\{-i \omega(k) t\}

Render the Waves

To render the waves, we can set the height of the vertex from the result of IFFT, calculated at each frame. This will be done easily. But we also want to calculate the normal of the vertex for proper lighting. This is done analytically by calculating the gradient of the height field:

\nabla h(\mathbf{x}, t) = \sum_{\mathbf{k}} i\mathbf{k} \tilde{h}(\mathbf{k}, t) \exp(i\mathbf{k} \cdot \mathbf{x})

Apply (I)FFT on GPU for Fast Computation

Suppose we have a grid with size 1024*1024, then we have to do 1024*1024 operations to get the height just for only 1 grid point! We have 1024*1024 grid points. Then the number of calculations would be 1024*1024*1024*1024! This is unacceptable. Fast fourier transform (FFT) is an efficient algorithm to compute the discrete Fourier transform (DFT) and its inverse. There're lot of tutorials explaining the algorithm itself. My short summary is that it utilizes some special properties of the terms in the summation, and apply Divide and Conquer to seperate the summation into even components and odd components, and so on and so forth. This reduce the complexity from $O(N^2)$ to $O(N \log N)$ .

There're also many variants of FFT algorithm. The one I used is called Stockham FFT, which is more suitable for parallel computation on GPU.

Figure 4. Comparison between Naive FFT and Stockham FFT. Stockham FFT has better memory access pattern for parallel computation on GPU. Instead of doing a bit reverse in the end like naive FFT, the pattern of indexing in each pass is the same in Stockham FFT.

In practice, I tabulated the twiddle factor and the indices of which 2 data to be combined, at the beginning of the FFT, which is a log(N) * N texture. So later in the FFT process, we can query the texture to get which two value we want to addup and what are their twiddle factor. The codes looks like this:

/// A sheet storing the idx and the complex exponential for stockham FFT. 
/// The sheet's shape is log(N) * N. Where N is the size of FFT. 
/// rg - twiddle factor. ba - indices of the two numbers to be combined.
[numthreads(1, 8, 1)]
void GenFFTParameters(uint3 id : SV_DispatchThreadID)
{
    // id.x is the current step(stage) of the FFT. 
    uint stage = id.x; 
    // It is calculating the distance between two numbers in the FFT. 
    uint distance = Size >> (stage + 1); 
    uint i = (2 * distance * (id.y / distance) + id.y % distance) % Size;
    float w = 2.0 * PI / (float) Size;
    float kn = (id.y / distance) * distance;
    float2 twiddle = C_Exp(float2(0, -w * kn));
    FFTParameterBuffer[id.xy] = float4(twiddle.x, twiddle.y, i, i + distance);
    FFTParameterBuffer[uint2(id.x, id.y + Size / 2)] = float4(-twiddle.x, -twiddle.y, i, i + distance);
}

For more information about Stockham FFT, you can refer to the CMU coursenote equation (13)^[5] and Multiprocessor FFTs by Paul N Swarztrauber equation (3.2)^[6]

Performance

The implementation can run at 600 FPS in 4K resolution on a 4090 GPU. The triangle count is 34.5k and verts count is 21.7k. The resolution for the FFT is 1024*1024.

Future work (if I have time...)

Try other spectrum models.
Add choppy waves computed from the transformation Jacobian.
The ocean's material now is not very satisfying to me. Lacking sense of transparency and refraction. Considering adding SSS material.
Offshore shallow water gradual change and underwater caustics.
Physical interactions.
Underwater scene rendering by Monzon et al. ^[7]

Some Other References Helped Me a Lot

Github repo by gasgiant: https://github.com/gasgiant/FFT-Ocean. The implementation used wave cascades to create a more detailed wave texture. The compute shader is written in a clear way that helps me a lot to understand.
Youtube tutorial by GetIntoGameDev: https://www.youtube.com/watch?v=sVps_gqlrqQ. This is a great tutorial on what is compute shader.