Flowing GPU compute water

Question

I have a civil engineering background and perform hydraulic and hydrologic analysis' regularly. They sell degrees for that sort of thing, but it's really not rocket science. I recently got it in my mind to implement the entire hydrologic and hydraulic process for a terrain on the GPU. I learned compute shaders only recently so I'm currently stuck being better at engineering than designing parallel GPU workflows.

You can calculate the amount of water generated during a rainfall event using the formula:
Q (CF/S) = c * I (in/hr) * A (acres)

I'm having difficulty moving beyond calculating the "acreage" of even the first area.

Current implementation overview:

1. The terrain is a regular grid of vertices at 1-unit intervals
2. A heightmap contains one R32 height value for each vertex
3. Currently, I'm only allowing flow in the 4 cardinal directions (no diagonals)
4. I'm using a Texture2D[int] as a stencil for vertices I have already analyzed

Current algorithm:

1. When the terrain tool has been active and is, now, not....
2. Clear the "stencil".
3. Scan the entire terrain for the lowest elevation.
4. That single point is the initial input to CS_Flood.
5. CS_Flood makes an X-axis pass.
6. Each input vertex is projected in both the X- and X+ directions up to 2048 times.
7. Finding an adjacent vertex with an OOB coordinate indicates the edge of the terrain in this direction. CurrentPoint is appended to the BoundaryPoints buffer and the projection loop for that direction is terminated. This was easy and works great every time.
8. Adjacent vertices with heights >= the current vertex's height are marked in the stencil and added to the NextPass buffer.
9. Adjacent vertices with heights < the current vertex's height indicates the peak of a ridge and terminates the projection loop. A future iteration of the flood fill might flow around the base of the ridge, up the "back" side of it and detect the same ridge a second time.
10. Any peak/ridge points that are detected more than once will not be a BoundaryPoint, for this purpose.
11. Any peak/ridge points detected exactly once are appended to BoundaryPoints and the projection loop in that direction is terminated.
12. CS_Flood makes a Z-axis pass with the same code, using the points generated by the X-axis pass as input.
13. Right now, CS_Flood continues alternating between the two directions indefinitely. Eventually, I'll terminate the overall loop whenever CS_Flood completes and the NextPass buffer is empty.

Ideally, at that point, BoundaryPoints would contain each vertex that occurs on the natural drainage divide. Water drops landing within the boundary eventually flow to the same low spot. Water drops landing ouside the boundary go "somewhere else".

Then:

1. Without clearing the stencil, re-scan the terrain for the lowest, non-stenciled, vertex.
2. Iterate CS_Flood.
3. Repeat until the stencil is full (or something similar).

3D is hard to perceive with these colors; this shows contour lines at integral elevations:
(a hole surrounded by a berm near the edge)

There are about 10 unique ways to drain across a vertex; giving each one a unique color looks like:
(visible circular tool marks, "ridges" show up nicely)

This shows every point generated by CS_Flood, boundary or otherwise, as a POINTLIST:

The algorithm always almost works. Sometimes, it even works correctly. Other times, the algorithm is clearly contained to the correct shape but will continue to output points indefinitely. As seen in the 3rd screenshot, sometimes it gets confused. There must be another situation/factor that I've overlooked. I would appreciate any help finding my oversight or suggestions of simpler and/or more elegant ways to attack the problem.

MissingPoint! can be included by band-aiding the algorithm to add every new BoundaryPoint detected to the NextPass buffer. During the next pass, 99% of the points generated by that band-aid will waste a small amount of GPU time determining that they can't go anywhere and doing nothing. During the first pass, sending LowestPoint along with the other NextPass points would also handle this specific scenario.

I know it's plausible and, given enough time, I'll be able to band-aid it enough to do what I want. I'd like to do it in a better, smarter, faster way, if possible, and I don't have enough experience yet to know better.

Answer 1

When a drop "tried" to visit a vertex, the stencil was marked with InterlockedExchange using the "original value" to determine if it was already stenciled (even though I just overwrote it).

The best algorithm I came up with gave the flood a "scratch pad" and a single rule: "don't flow down-hill" (heights equal to or greater). That eliminated almost all of the complicated tests. While it is generally good at not flowing over peaks/ridges, it does flow along them because the adjacent vertices are "flat". This occasionally allows drops to sneak past ridge-lines.

Each of the "too far" points is, then, "flowed" and will flow "into" the drainage area (stops at a 1) or not (stops at a 0). The "nots" are discarded and the corrected scratch-pad is copied into the "final". If the final is already stenciled, the scratch pad is discarded. (Future: those collisions should collectively represent the outer boundary of the current drainage area.)

At 10FPS:

The "nots" are shown in red, once the large area is copied into the final and becomes green, then the algorithm repeats for the remaining unstenciled areas.