Strange vibrational ripples in my shallow water

Question

I tried to implement shallow water equations in Unity, but I ran into a strange error. I get these strange ripples in my water. I took some screenshots:

And the video you can find here: https://www.youtube.com/watch?v=crXLrvETdjA

I based my code on paper. Quick modeling and visualization of hydraulic erosion on the Xing Mei GPU . And here you can find the complete solution code: http://pastebin.com/JktpizHW (or see below). Each time I use the formula from the article, I added my number as a comment.

I tried different timestamps, for the video I used 0.02, its decrease simply made it oscillate more slowly. I also tried a large grid (the video uses 100, I tried 200, but then the ripples were smaller.) I checked all the formulas several times and cannot find any errors.

Anyone who can understand what is going wrong?

Additional Information:

As you can see from pastebin, I programmed it in C #. I used Unity as my render engine, and I just use a grid mesh to render water. I am changing the vertex y element of the grid according to the calculated height.

The DoUpdate method gets the parameter float[][] lowerLayersHeight , which is basically the height of the terrain under water. In the video, it's just 0 .

 public override void DoUpdate(float dt, float dx, float[][] lowerLayersHeight) { int x, y; float totalHeight, dhL, dhR, dhT, dhB; float dt_A_g_l = dt * _A * g / dx; //all constants for equation 2 float K; // scaling factor for the outflow flux float dV; for (x=1 ; x <= N ; x++ ) { for (y=1 ; y <= N ; y++ ) { // // 3.2.1 Outflow Flux Computation // -------------------------------------------------------------- totalHeight = lowerLayersHeight[x][y] + _height[x][y]; dhL = totalHeight - lowerLayersHeight[x-1][y] - _height[x-1][y]; //(3) dhR = totalHeight - lowerLayersHeight[x+1][y] - _height[x+1][y]; dhT = totalHeight - lowerLayersHeight[x][y+1] - _height[x][y+1]; dhB = totalHeight - lowerLayersHeight[x][y-1] - _height[x][y-1]; _tempFlux[x][y].left = Mathf.Max(0.0f, _flux[x][y].left + dt_A_g_l * dhL ); //(2) _tempFlux[x][y].right = Mathf.Max(0.0f, _flux[x][y].right + dt_A_g_l * dhR ); _tempFlux[x][y].top = Mathf.Max(0.0f, _flux[x][y].top + dt_A_g_l * dhT ); _tempFlux[x][y].bottom = Mathf.Max(0.0f, _flux[x][y].bottom + dt_A_g_l * dhB ); float totalFlux = _tempFlux[x][y].left + _tempFlux[x][y].right + _tempFlux[x][y].top + _tempFlux[x][y].bottom; if (totalFlux > 0) { K = Mathf.Min(1.0f, _height[x][y] * dx * dx / totalFlux / dt); //(4) _tempFlux[x][y].left = K * _tempFlux[x][y].left; //(5) _tempFlux[x][y].right = K * _tempFlux[x][y].right; _tempFlux[x][y].top = K * _tempFlux[x][y].top; _tempFlux[x][y].bottom = K * _tempFlux[x][y].bottom; } //swap temp and the real one after the for-loops // // 3.2.2 Water Surface // ---------------------------------------------------------------------------------------- dV = dt * ( //sum in _tempFlux[x-1][y].right + _tempFlux[x][y-1].top + _tempFlux[x+1][y].left + _tempFlux[x][y+1].bottom //minus sum out - _tempFlux[x][y].right - _tempFlux[x][y].top - _tempFlux[x][y].left - _tempFlux[x][y].bottom ); //(6) _tempHeight[x][y] = _height[x][y] + dV / (dx*dx); //(7) //swap temp and the real one after the for-loops } } Helpers.Swap(ref _tempFlux, ref _flux); Helpers.Swap(ref _tempHeight, ref _height); } 

Answer 1

I fixed it myself! Although about this while driving to a friend. The problem is quite simple, what I do in the code I listen to, for each cell (or grid point) it calculates the flow, then the height, and then I move on to the next cell. I must first calculate the flow for all cells, then repeat a second time across all cells and calculate their height. Thus, the code becomes:

 for (x=1 ; x <= N ; x++ ) { for (y=1 ; y <= N ; y++ ) { // // 3.2.1 Outflow Flux Computation // -------------------------------------------------------------- *** } } for (x=1 ; x <= N ; x++ ) { for (y=1 ; y <= N ; y++ ) { // // 3.2.2 Water Surface // --------------------------------------------------------------------------- *** } } Helpers.Swap(ref _tempFlux, ref _flux); Helpers.Swap(ref _tempHeight, ref _height); 

(Of course *** becomes the corresponding code from the question above).

Now I have a simulation of working water.