Waves, finally

Before I ended up with a cube floating stably on perfectly flat water at $y=0$. Time to give it something to actually float on. This post has two halves, build a wave system in C# that both the visualizer and the buoyancy code agree on, and rewire the DLL so it stops pretending the water is flat.

The architectural pivot

I originally planned for the C++ DLL to know about waves itself, take a WaveParams struct over the P/Invoke boundary, evaluate $W(x,z,t)$ in C++, classify triangles, etc. After actually writing some of it I changed my mind. The DLL has no business knowing what kind of waves we're running. Its job is "given a triangle and the water heights at each corner, give me hydrostatic force and torque." That's a pure geometric primitive, wave shape doesn't enter into it. This makes it less coupled and more flexible for other people to use.

So the new split is C# samples wave heights at every mesh vertex's world (x,z) position and passes a float[] of heights into the DLL alongside the transform matrix. The DLL looks at the heights, classifies triangles, clips against the per-vertex surface, and returns the same six floats it always did.

Reasons this is better:

1. Swap wave models without recompiling the DLL. We decouple the wave from the DLL.
2. The clipper already does linear interpolation for clip-point positions. Lerping the water heights between two endpoint heights uses the same parameter $t$. Zero new approximation introduced by passing data instead of evaluating the function.
3. Marshalling cost is trivial. ~24 cube vertices × 4 bytes = 96 bytes per frame. Not measurable.
4. Future Burst job parallelises the height sampling across cores for free, the single-threaded DLL doesn't get that. Could be done with OpenMP or similar but takes more effort.

Wave function, Catlike's Gerstner

Could've written my own, but Jasper Flick at Catlike Coding already has a Gerstner waves tutorial that's clean, correct, and licensed MIT-0. So I copied it. The C# and shader implement the same Gerstner function so the physics inversion always works on the wave that's actually visible. This feels like duplication but in reality we are forced to because the GPU handles the visuals and the CPU the physics. We could have a compute shader run the calculations but that would lead to a head ache trying to sync everything together so figured that for this I would just split it.

A single Gerstner wave displaces a rest position $\mathbf{p}_0$ to a new world point:

$\mathbf{p} = \mathbf{p}_0 + \begin{bmatrix} D_x \cdot a \cos f \\ a \sin f \\ D_y \cdot a \cos f \end{bmatrix}, \quad f = k(\mathbf{D} \cdot \mathbf{p}_{0,xz} - c t), \quad a = \frac{s}{k}$

with $k = 2\pi/\lambda$, $c = \sqrt{g/k}$, $\mathbf{D}$ the unit horizontal direction, and steepness $s \in [0,1)$. Sum three of them with different directions/wavelengths and you get a passable ocean.

C# port:

static Vector3 GerstnerDisp(Vector4 wave, float px, float pz, float t)
{
    float steepness = wave.z;
    float wavelength = wave.w;
    float k = 2f * Mathf.PI / wavelength;
    float c = Mathf.Sqrt(G / k);
    Vector2 d = new Vector2(wave.x, wave.y).normalized;
    float f = k * (d.x * px + d.y * pz - c * t);
    float a = steepness / k;
    return new Vector3(d.x * (a * Mathf.Cos(f)),
                               a * Mathf.Sin(f),
                       d.y * (a * Mathf.Cos(f)));
}

The Gerstner sampling problem

There's a subtle catch I had to deal with. Gerstner displaces in $x$ and $z$ as well as $y$. So GerstnerDisp(p_0) returns the displacement of the surface particle that originated at $\mathbf{p}_0$. For buoyancy I don't want that, I want the height of the surface above a given world $(x, z)$. They're different points.

It's an inverse mapping problem. Given world $(x, z)$, find the rest position $\mathbf{p}_0$ such that

$\mathbf{p}_0 + \Delta_{xz}(\mathbf{p}_0) = (x, z)$

then evaluate the height at that rest position. Closed-form? Doesn't exist. Fixed-point iteration? Works great, converges in 3-4 iterations for any reasonable steepness:

public float SampleHeight(float x, float z, float t)
{
    float px = x, pz = z;
    for (int i = 0; i < sampleIterations; i++)
    {
        Vector3 d = GerstnerDisp(waveA, px, pz, t)
                  + GerstnerDisp(waveB, px, pz, t)
                  + GerstnerDisp(waveC, px, pz, t);
        px = x - d.x;
        pz = z - d.z;
    }
    Vector3 disp = GerstnerDisp(waveA, px, pz, t)
                 + GerstnerDisp(waveB, px, pz, t)
                 + GerstnerDisp(waveC, px, pz, t);
    return disp.y;
}

For my default steepness of 0.25 and 60m wavelength the horizontal drift is small enough that even one iteration would be visually fine. But running four costs nothing per cube and means the per-vertex heights match the visible surface to floating-point precision. The water mesh and the boats line up. This however I realized later wasn't the case for more complex meshes.

Wave manager

WaveManager is a MonoBehaviour singleton owning the three Gerstner waves as Vector4 fields. It does two things every frame, pushes the wave parameters into the water material (so the shader uses identical math), and exposes SampleHeight(x, z, t) for the buoyancy controller to call per vertex.

void OnValidate() { if (waterMaterial != null) PushToMaterial(); }

void PushToMaterial()
{
    waterMaterial.SetVector(WaveAId, waveA);
    waterMaterial.SetVector(WaveBId, waveB);
    waterMaterial.SetVector(WaveCId, waveC);
}

Critical that there's exactly one source of truth here. If the shader and the C# code use different parameters for one frame, the cube floats above or below the visible surface in a really obvious way.

URP shader, see-through with depth fade

Catlike's tutorial ships a built-in-pipeline surface shader. My project is URP, so the shader needed porting. The math (GerstnerWave HLSL function, vertex displacement, tangent/binormal accumulation, normal cross product) is line-for-line identical to Catlike's. Only the boilerplate changes, CGPROGRAM to HLSLPROGRAM, appdata_full to custom Attributes, surf function to fragment shader calling UniversalFragmentPBR.

One important fix though, Catlike's shader computes waves in object space. That works fine for a stationary plane at the origin, but if you ever rotate or move the water, the waves slide with it, and the buoyancy code (which samples in world space) desyncs. Fix is one line:

float3 gridPointWS = TransformObjectToWorld(IN.positionOS.xyz);
// ... waves computed from gridPointWS, not positionOS

While I had the shader open I added two URP niceties since it's worthless for my demo if what's happening underwater isn't visible:

• Transparency, Tags { "Queue"="Transparent" }, Blend SrcAlpha OneMinusSrcAlpha, ZWrite Off. Material gets an _Opacity slider.
• Depth fade, sample the opaque scene depth behind the water with SampleSceneDepth, blend toward fully opaque as the depth difference grows.

float2 uv = IN.screenPos.xy / IN.screenPos.w;
float sceneDepth = LinearEyeDepth(SampleSceneDepth(uv), _ZBufferParams);
float surfDepth  = IN.screenPos.w;
float depthDiff  = sceneDepth - surfDepth;
float depthAlpha = saturate((depthDiff - _DepthFadeNear) / (_DepthFadeFar - _DepthFadeNear));
float alpha = _Opacity + (1.0 - _Opacity) * depthAlpha;

Now I can see the cube sinking through the surface, settling at its waterline, half-submerged with the rest visible underwater. Much easier to debug than opaque blue.

Wiring heights into the DLL

Two changes in BuoyancyController.cs. First, allocate the height array and the matrix array once in Start() instead of every frame:

localVerts = verts;
vertexHeights = new float[verts.Length];
float[] matrix = new float[16];

Second, sample heights in FixedUpdate and pass them to the DLL:

float t = Time.fixedTime;
WaveManager wm = WaveManager.Instance;
for (int i = 0; i < localVerts.Length; i++)
{
    Vector3 wv = m.MultiplyPoint3x4(localVerts[i]);
    vertexHeights[i] = wm.SampleHeight(wv.x, wv.z, t);
}
NativeBridge.ComputeBuoyancy(handle, matrix, vertexHeights, vertexHeights.Length, res);

Inside the DLL, accumulateBuoyancy now takes per-corner heights (ha, hb, hc), averages them for the centroid water height, and oneAboveSplit / oneBelowSplit interpolate heights at the clip points using the same parameter $t$ they already used for positions:

float step1 = glm::clamp(dFirst / (dFirst - dAbove), 0.0f, 1.0f);
glm::vec3 p1 = glm::mix(first, above, step1);
float    hp1 = glm::mix(hFirst, hAbove, step1);

The clamp matters. Without it, when the wave crest pushes a vertex's water height up to where it nearly equals the vertex's $y$, the denominator goes near zero and step blows up to thousands. Clip points end up kilometers from the cube, sub-triangle areas explode, and you get a single-frame force of 210 MN that launches the cube into orbit, happend to me...

Things that went wrong

Forgot the DLL was wave-agnostic the first time, spent an hour wondering why turning waves on changed nothing for the buoyancy. Turns out the DLL wasn't actually using the height array I was passing, the helpers still had their literal 0.0f arguments inside ComputeBuoyancy's switch statement. Whole point of the architectural pivot defeated by one missed line.

Parameter-order bug between forward declaration and call site, refactored accumulateBuoyancy to take per-corner heights, changed the signature to (... rho, gravity, ha, hb, hc, ...). Forgot to update the forward declaration at the top of the file, which had (... ha, hb, hc, rho, gravity, ...). C++ happily compiled because the types match (all float). At runtime, rho got the value hb = 0, so the force computation became 0 * gravity * ... and the cube just sank straight down, no force, no torque, no warning. Caught it by adding a debug log that printed the cube's position, the first vertex's height, and the depth, and noticed depth was clearly negative but force was still 0. Lesson learned, when a function takes a lot of float parameters, swap two and the compiler doesn't care.

Cube flew across the map first time waves were enabled. This was the explosion described above, before the clamp went in. With waves at zero amplitude everything had been fine, turn steepness up to 0.25 for the waves and a single bad clip point per frame is enough to send momentum to infinity. Quick fix (clamp(...)) once I'd figured out which value was nonsense.

Wish I had taken a video before fixing this... Would have been great for the blog.

Result