I’ve been experimenting ab it for the past week and have found an open source C++ implementation of the Shallow Water equations. I’ve come across these equations quite a lot when searching for water simulation in games. I’ve managed to implement it in UE4 in C++, having changed a few things such as variable types and names to work with the engine. While it’s roughly working, the height values eventually build up and reaches into infinity very quickly. If I use them to drive the Z coordinates of vertices, the engine crashes with a mesh out of bounds error (of course caused by the numbers going towards infinity). Unfortunately, I’ve only scratched the surface of understanding equations, so trying to find a solution for it is way beyond my math skills.
The implementation is below:
void AWaterSurface::UpdateHeight()
{
float Temp;
int x;
int y;
x = 0;
// Set border cells
for (x = 0; x < SurfaceVertices + 2; x++)
{
H[IX(x, 0)] = H[IX(x, 1)];
U[IX(x, 0)] = U[IX(x, 1)];
V[IX(x, 0)] = -V[IX(x, 1)];
H[IX(x, SurfaceVertices + 1)] = H[IX(x, SurfaceVertices)];
U[IX(x, SurfaceVertices + 1)] = U[IX(x, SurfaceVertices)];
V[IX(x, SurfaceVertices + 1)] = -V[IX(x, SurfaceVertices)];
}
for (y = 0; y < SurfaceVertices + 2; y++)
{
H[IX(0, y)] = H[IX(1, y)];
U[IX(0, y)] = -U[IX(1, y)];
V[IX(0, y)] = V[IX(1, y)];
H[IX(SurfaceVertices + 1, y)] = H[IX(SurfaceVertices, y)];
U[IX(SurfaceVertices + 1, y)] = -U[IX(SurfaceVertices, y)];
V[IX(SurfaceVertices + 1, y)] = V[IX(SurfaceVertices, y)];
}
for (y = 0; y < SurfaceVertices; y++)
{
for (x = 0; x < SurfaceVertices + 1; x++)
{
Temp = (U[IX(x + 1, y + 1)] - U[IX(x, y + 1)]);
Hx[IX(x, y)] = (H[IX(x + 1, y + 1)] + H[IX(x, y + 1)]) / 2 - (Dt / (2 * Dx)) * Temp;
Temp = (U[IX(x + 1, y + 1)] * U[IX(x + 1, y + 1)] / H[IX(x + 1, y + 1)] + (G / 2) * H[IX(x + 1, y + 1)] * H[IX(x + 1, y + 1)] -
(U[IX(x, y + 1)] * U[IX(x, y + 1)] / H[IX(x, y + 1)] + G / 2 * H[IX(x, y + 1)] * H[IX(x, y + 1)]));
Ux[IX(x, y)] = (U[IX(x + 1, y + 1)] + U[IX(x, y + 1)]) / 2 - (Dt / (2 * Dx)) * Temp;
Temp = (U[IX(x + 1, y + 1)] * V[IX(x + 1, y + 1)] / H[IX(x + 1, y + 1)] -
U[IX(x, y + 1)] * V[IX(x, y + 1)] / H[IX(x, y + 1)]);
Vx[IX(x, y)] = (V[IX(x + 1, y + 1)] + V[IX(x, y + 1)]) / 2 - (Dt / (2 * Dx)) * Temp;
}
}
for (y = 0; y < SurfaceVertices + 1; y++)
{
for (x = 0; x < SurfaceVertices; x++)
{
Temp = (V[IX(x + 1, y + 1)] - V[IX(x + 1, y)]);
Hy[IX(x, y)] = (H[IX(x + 1, y + 1)] + H[IX(x + 1, y)]) / 2 - (Dt / (2 * Dy)) * Temp;
Temp = (V[IX(x + 1, y + 1)] * U[IX(x + 1, y + 1)] / H[IX(x + 1, y + 1)] -
V[IX(x + 1, y)] * U[IX(x + 1, y)] / H[IX(x + 1, y)]);
Uy[IX(x, y)] = (U[IX(x + 1, y + 1)] + U[IX(x + 1, y)]) / 2 - (Dt / (2 * Dy)) * Temp;
Temp = (V[IX(x + 1, y + 1)] * V[IX(x + 1, y + 1)] / H[IX(x + 1, y + 1)] + (G / 2) * H[IX(x + 1, y + 1)] * H[IX(x + 1, y + 1)] -
(V[IX(x + 1, y)] * V[IX(x + 1, y)] / H[IX(x + 1, y)] + (G / 2) * H[IX(x + 1, y)] * H[IX(x + 1, y)]));
Vy[IX(x, y)] = (V[IX(x + 1, y + 1)] + V[IX(x + 1, y)]) / 2 - (Dt / (2 * Dy)) * Temp;
}
}
for (y = 1; y < SurfaceVertices + 1; y++)
{
for (x = 1; x < SurfaceVertices + 1; x++)
{
Temp = (Dt / Dx) * (Ux[IX(x, y - 1)] - Ux[IX(x - 1, y - 1)]) + Dt / Dy * (Vy[IX(x - 1, y)] - Vy[IX(x - 1, y - 1)]);
H[IX(x, y)] = H[IX(x, y)] - Damping * Temp;
U[IX(x, y)] = U[IX(x, y)] - (Dt / Dx) * (Ux[IX(x, y - 1)] * Ux[IX(x, y - 1)] / Hx[IX(x, y - 1)] + (G / 2) * Hx[IX(x, y - 1)] * Hx[IX(x, y - 1)] -
(Ux[IX(x - 1, y - 1)] * Ux[IX(x - 1, y - 1)] / Hx[IX(x - 1, y - 1)] + (G / 2) * Hx[IX(x - 1, y - 1)] * Hx[IX(x - 1, y - 1)]))
- (Dt / Dy) * (Vy[IX(x - 1, y)] * Uy[IX(x - 1, y)] / Hy[IX(x - 1, y)] -
(Vy[IX(x - 1, y - 1)] * Uy[IX(x - 1, y - 1)] / Hy[IX(x - 1, y - 1)]));
V[IX(x, y)] = V[IX(x, y)] - (Dt / Dx) * (Ux[IX(x, y - 1)] * Vx[IX(x, y - 1)] / Hx[IX(x, y - 1)] -
(Ux[IX(x - 1, y - 1)] * Vx[IX(x - 1, y - 1)] / Hx[IX(x - 1, y - 1)]))
- (Dt / Dy) * (Vy[IX(x - 1, y)] * Vy[IX(x - 1, y)] / Hy[IX(x - 1, y)] + (G / 2) * Hy[IX(x - 1, y)] * Hy[IX(x - 1, y)] -
(Vy[IX(x - 1, y - 1)] * Vy[IX(x - 1, y - 1)] / Hy[IX(x - 1, y - 1)] + (G/2) * Hy[IX(x - 1, y - 1)] * Hy[IX(x - 1, y - 1)]));
}
}
}
As far as I know, “H” is water height, and “U” and “V” refer to the flow velocity in the X and Y directions. I eventually want these three values to go into a render target to make use of the GPU to allow for larger water surfaces, and to use the flow velocity to control the direction and speed of the ripple normal maps.
EDIT: the IX(x,y) is a macro that allows me to search a 1D array as if it’s a 2D array.