Dziwne oscylujące fale w mojej płytkiej implementacji

Question

Próbowałem wprowadzić równania płytkiej wody w Unity, ale uruchomiłem dziwny błąd. Dostaję te dziwne oscylujące fale w mojej wodzie. Zrobiłem jakiś zrzut ekranu:

i wideo można znaleźć tutaj: https://www.youtube.com/watch?v=crXLrvETdjA

oparłem mój kod na papierze Fast Hydraulic Erosion Simulation and Visualization on GPU przez Xing Mei. Możesz znaleźć cały kod solver tutaj: http://pastebin.com/JktpizHW (lub patrz poniżej). Za każdym razem, gdy używam formuły z artykułu, dodałem jej numer jako komentarz.

Próbowałem różnych timesteps, dla wideo użyłem 0.02, obniżając go po prostu spowodował, że będzie wolniej oscylować. Próbowałem także większej siatki (wideo używa 100, próbowałem 200, ale potem zmarszczki były po prostu mniejsze.) Sprawdziłem wszystkie formuły kilka razy i nie mogłem znaleźć żadnego błędu.

Ktoś tutaj może dowiedzieć się, co się dzieje?

Dodatkowe info:

Jak widać z pastebin, ja go programować w C#. Użyłem Unity jako mojego silnika do wizualizacji i używam tylko siatki do wizualizacji wody. Zmieniam składnik wierzchołka siatki y, dopasowany do wysokości, którą obliczam.

Metoda DoUpdate otrzymuje parametr float[][] lowerLayersHeight, który jest w zasadzie wysokością terenu pod wodą. Na filmie jest tylko 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

Naprawiłem to sam! Chociaż to podczas jazdy do znajomego. Problem jest dość prosty, to, co robię w błędnym kodzie, jest dla każdej komórki (lub punktu siatki) obliczyć strumień, następnie wysokość, a następnie przejść do następnej komórki. Powinienem najpierw obliczyć strumień dla wszystkich komórek, następnie powtórzyć po raz drugi we wszystkich komórkach i obliczyć ich wysokość. Tak więc staje się kod:

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); 

(Oczywiście *** staje się odpowiedni kod z pytania powyżej.)

Teraz mam symulacji wody roboczej.