Limited water simulation

Development Rendering
1

For a serious project we’re planning to work on with a team, we need to have a limited form of water simulation. Right now we’re exploring a lot of options, but water is of course extremely complex and we are not exceptionally skilled at fluid dynamics. However, we feel that we can make at least a basic form of simulation. We are exploring NVIDIA Flex and Cataclysm, but those are of course still experimental and we would have to compile the engine source, so we are not sure it is suitable for production work, at least not until the latest DirectCompute versions are released.

We do not need a full-blown water simulation like the aforementioned NVIDIA systems. What we need at minimum is a way of a water surface rising to fill an area based on the height of a “source” or ocean and then more-or-less realistically propagating to fill the whole area. When a wall is removed, the water surface should propagate to fill the next area. We believe the system used in Cities Skylines is essentially similar.

Right now, we’re thinking of dividing the world into cells, with each cell effectively holding the water height and the flow direction in the X and Y directions. These cells would then be used to modify the Z positions of vertices of a procedural plane mesh. We then need to detect enclosed areas to avoid filling isolated areas (as in those protected by dikes, dams and other barriers). In principle we would use the flow direction of a cell to create a dynamic flow map that adjusts the direction of a normal map to fake the flow of water. The amount of cells determine the resolution and complexity of the simulation. We’re also thinking of exploring cellular automata to attempt to simulate the flow.

So what I’m asking here essentially if this is a reasonable direction to go and if anyone more experienced with this has some pointers? We have done extensive research across the internet for the past few days, but almost all attempts result in articles about a full fluid simulation or a 3D/2D grid-based simulation (like Minecraft or Dwarf Fortress). Other results that come close have little to no explanation on how to implement something similar.

2

For me, everything you mentioned is just a simple water, with standard normal texture panner flow. There is no need to make some specific directions and additional effects, if it not affects gameplay.
So, use simple panner and spline for that.

3

I think you misunderstand me. Look at how Cities Skylines is handling water. A simple flat surface is not enough, the water has to appear to flood an area that’s below water level. Much like dynamic ripples created by the old FluidSurfaceActor and what Epic has demonstrated with Blueprint Render Targets, but that it should also raise the water level directly, as in the vertices moving up. And it does affect gameplay, that’s the point. We have to simulate at least somewhat realistically a flooding if a dam breaks (hence the Skylines reference)

4

A dynamic ripples shader and a couple RenderTargets could do what you need. Combine it with a heightmap of the flow limiting geometry using a SceneCapture2D. Use orthographic projection and apply a Post Process material that returns scene depth. If you don’t need solid collision, not much point in doing it CPU side and using a procedural mesh component. It would be dreadfully slow. Using Render Target Lookup on a per interacting object basis would be faster.

5

What’s the level of interaction with the water or is this purely an aesthetic thing? For instance if the dam breaks does the city floods the exact same way every time with some minor building interactions?

6

It needs to be dynamic, though not extremely detailed or realistic. This is because things might be different, such as the landscape, location of the breach and placement of buildings. Although the water doesn’t necessarily have to interact with the buildings, it doesn’t have to be that accurate.

7

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.

8

Usually there’s a dispersion factor like 0.99 or 0.999, Take this example on shadertoy. Look in buffer B for the line: nu = 0.99*nu; I noticed your code is lacking that.

What sort of performance are you getting? I still imagine it would be faster to do GPU side and read the pixels.

9

The performance is quite good, but I’m testing with a relatively small and low-res plane. I do intend to send the values to a texture later on when it works.

If I understand correctly, as I’m still wrapping my head around it, dispersion results in the waves dissipating over time, so they eventually disappear. In my code, there’s already a damping factor. If I decrease that value, the waves do dissipate (disperse, am I correct?) much faster. However, then when I “add” water by increasing the water level of cells, the infinity problem also happens much faster. I double-checked my code with the original source, and there are no typos. It appears to be caused by the border cells, as they reflect the waves back, but then they won’t stop. I’ll keep experimenting a little with the border cells to see if that could be the problem.

This code selects a random cell and adds a value to it (set in the editor).



void AWaterSurface::CreateSource()
{
	int DesiredX = (int)(FMath::FRandRange(1, SurfaceVertices + 1));
	int DesiredY = (int)(FMath::FRandRange(1, SurfaceVertices + 1));

	H[IX(DesiredX, DesiredY)] += SourceStrength;
}