Strange oscillating ripples in my shallow water implementation <p>I've been trying to implement the shallow water equations in Unity, but I've run in a weird bug. I get these strange oscillating ripples in my water. I made some screenshot:</p> <p><img src="https://i.stack.imgur.com/9FcuF.jpg" alt="enter image description here"> <img src="https://i.stack.imgur.com/JpCOk.jpg" alt="enter image description here"></p> <p>And a video you can find here: <a href="https://www.youtube.com/watch?v=crXLrvETdjA" rel="nofollow noreferrer">https://www.youtube.com/watch?v=crXLrvETdjA</a></p> <p>I based my code on the paper <a href="http://evasion.imag.fr/Publications/2007/MDH07/FastErosion_PG07.pdf" rel="nofollow noreferrer">Fast Hydraulic Erosion Simulation and Visualization on GPU</a> by Xing Mei. And you can find the entire solver code here: <a href="http://pastebin.com/JktpizHW" rel="nofollow noreferrer">http://pastebin.com/JktpizHW</a> (or see below.) Each time I use a formula from the paper, I added its number as a comment.</p> <p>I tried different timesteps, for the video I used 0.02, lowering it just made it oscillate slower. I also tried a bigger grid (video uses 100, I tried 200 but then the ripples just were smaller.) I checked all formulas several times, and can't find any error.</p> <p>Anyone here who can figure out what's going wrong?</p> <p>Extra info:</p> <p>As you can see from the pastebin, I programmed it in c#. I used Unity as my engine for visualisation, and I'm just using a grid mesh to visualise the water. I alter the mesh's vertex y component to match the height I calculate.</p> <p>The DoUpdate method gets a <code>float[][] lowerLayersHeight</code> parameter, that's basically the height of the terrain under the water. In the video it's just all <code>0</code>.</p> <pre><code>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 outﬂow ﬂux float dV; for (x=1 ; x &lt;= N ; x++ ) { for (y=1 ; y &lt;= N ; y++ ) { // // 3.2.1 Outﬂow 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 &gt; 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); } </code></pre>