By popular demand, here's the follow up to my post about the water tech in my game, Days on the Lake. You can read the first part here: https://ldjam.com/events/ludum-dare/53/days-on-the-lake/simulating-water-in-days-on-the-lake

I'm going to start by showing y'all the simulation code (in HLSL), and then I'll go section by section and explain what's happening. Note: as written this can't be copy-pasted into an HLSL shader. HLSL is C-like and requires forward declaration of functions, and I've reordered a lot of the code here to make it more readable. As well, I don't show my noise implementation. I recommend reading about simple GPU noise implementations here: https://thebookofshaders.com/11/
```
float2 frag(v2fcustomrendertexture IN) {
// Sample textures
float4 terrainHeight = tex2D(TerrainTex, IN.localTexcoord);
float4 currentWater = tex2D(_SelfTexture2D, uv);
// Calculate derivatives and wave speed
float4 d = getDerivatives(IN.localTexcoord);
float2 cSquared = getWaterC(derivatives);
// Calculate forces
float forcing = getForcing(IN.localTexcoord);
float restore = -currentWater.x * _Restore;
float damping = -currentWater.y * _Damping;
float wave = cSquared.x * d.ddx2 + cSquared.y * d.ddy2 + (d.ddx * d.ddx + d.ddy * d.ddy) * _FakeFactor;
// Calculate next wave state
float terrainHeightClamp = _TerrainHeightOffset - terrainHeight.r * _TerrainHeightScale;
float nextVelocity = currentWater.y + (wave + forcing + restore + damping) * _DeltaTime;
float nextHeight = clamp(currentWater.x + v * _DeltaTime, -terrainHeightClamp, TerrainHeightClamp);
return float4(nextHeight, nextVelocity, 0, 1);
}
Derivatives getDerivatives(float2 uv) {
Derivatives derivatives;
float texelX = 1 / _CustomRenderTextureWidth;
float texelY = 1 / _CustomRenderTextureHeight;
float4 xy1 = tex2D(_SelfTexture2D, uv);
float4 x0 = tex2D(_SelfTexture2D, uv + float2(-texelX, 0));
float4 x2 = tex2D(_SelfTexture2D, uv + float2(texelX, 0));
derivatives.ddx = (x2.x - x0.x) / (2 * _GridSize);
derivatives.ddx2 = (x0.x + -2 * xy1.x + x2.x) / (_GridSize * _GridSize);
float4 y0 = tex2D(_SelfTexture2D, uv + float2(0, -texelY));
float4 y2 = tex2D(_SelfTexture2D, uv + float2(0, texelY));
derivatives.ddy = (y2.x - y0.x) / (2 * _GridSize);
derivatives.ddy2 = (y0.x + -2 * xy1.x + y2.x) / (_GridSize * GridSize);
return derivatives;
}
float2 getWaterC(Derivatives d) {
float wavelengthX = 5 * sqrt(abs(d.ddx2)) + abs(d.ddx);
float wavelengthY = 5 * sqrt(abs(d.ddy2)) + abs(d.ddy);
float cxSquared = (_Gravity * wavelengthX) / (2 * PI);
float cySquared = (_Gravity * wavelengthY) / (2 * PI);
return float2(cxSquared, cySquared);
}
float getForcing(float2 uv) {
float forcing = valueNoise(float3(uv.x * _ForcingFreq, uv.y * _ForcingFreq, _Time.z * _ForcingTimescale)) - .5;
float s = sign(forcing);
forcing = forcing * forcing * s;
return forcing * _Forcing;
}
struct Derivatives {
float ddx;
float ddx2;
float ddy;
float ddy2;
}
```
This is a lot! To we'll look at frag for the basic outline of the method. First, we sample the terrain height and the last water texture at the current pixel. Then we get the derivatives (the du/dx stuff) of the water texture, and calculate wave speed with them. Then we calculate all of the "forces" acting on our water. These are "forces" in a loose sense, but it's a reasonable way to think about them. Then we integrate forward with time, and return it as the next texture step!
In this, I'm integrating forward in time with the forward euler method. The forward euler method is inaccurate and doesn't conserve energy, but it's sufficient for a simple simulation in a game. Simulations with this method take the form you see in the code: nextState = currentState + currentRateOfChange * deltaTime. There are more accurate things you can do here. The midpoint method in time is often better, or a variation of verlet integration. You could also do one of the runge-kutta family of methods, or even an iterative method like the backward euler method. I won't go into these here, but they're good things to read up on if you're interested in simulating things!
In the getDerivatives function, I use the central difference method to get derivatives of the water height field in space. The central difference method is the simplest method to do this. I don't use the built in ddx and ddy functions here because I can't construct second derivatives with them, because of technical limitations. In order to make consistent first and second derivatives, I just manually calculate both. The built in ddx and ddy take advantage of the GPU's parallel execution over many pixels to calculate derivatives quickly, but they're only able to look at 2x2 blocks. This is insufficient to calculate a second derivative, sadly!
Moving on, we have the getWaterC function. This estimates the wave speed of water given the first and second derivatives. I basically played around in graphing software to figure out something that estimated the wavelength well, because we need it to calculate C. What I landed on wasn't great, but was sufficient. It varies a lot over different parts of a wave, but kind of hovers around the actual wavelength. You can read up on the calculation of the actual wave speed, given a wavelength here: https://en.wikipedia.org/wiki/Dispersion(waterwaves). I probably should have implemented the shallow water effects as well at least. It'll be an interesting thing to investigate going forward, as well as the nonlinear large wave effects in deep water.
Next, we have the "forces". I keep putting "forces" in quotes because they're not really forces as you see in rigidbody simulation. I'm not considering the effect of these in relation to a mass, I'm just using them as terms to derive an acceleration. The term "force" is good for intuition though, so I like it. The "forces" I use are damping, restore, forcing, and wave.
damping is just this way of bleeding off velocity, to what would be friction. It's used in a lot of simulations and is usually taken as proportional to velocity.
The restore term is the way I force water to get back to the "normal" water height, and I take it as proportional to height. This is a fictitious force, but I find it helpful in stabilizing these simulations. You can imagine it's a spring pulling each pixel back to center.
The forcing term is what I use to generate the idle motion of water. I won't go too deep into the getForcing function, but it basically samples some value noise and squares it to make the peaks sharper. This is a place you can really get creative! Copy the code and try weird things!
The wave term is the meat of the simulation. It's what makes the waves propagate like real waves! It can be understood in two parts. First is the proper simple wave equation, d^2u/dt^2 = c^2 * d^2u/dx^2. This just says the amount the wave wants to accelerate up or down is proportional to the second derivative of it in space, or how much it's smiling or frowning. Frowning waves want to go down, smiling waves want to go up. The second part is the part I added (dx * dx + dx * dy) * _FakeFactor. This, as the factor's name implies, is fake, and just me trying to make the waves look more interesting. It basically says steep water waves want to stay "up". It's overpowered by the proper wave equation part, but I like to think it does something.
After that, we just integrate forward in time as I've talked about, and return. The clamping to the shore part is something I added later, to make the edges look better.
This has gotten pretty long, so I think I'll end this one here. I have one more post to do about the way I shade the water, so I hope to see folks back for that!