Fluid

/// <summary>
/// Update the fluid. Part 3.2 in the paper.
/// </summary>
/// <param name="dt">delta time</param>
/// <param name="dx">delta x. Distance between middles of grid-cells.</param>
/// <param name="lowerLayersHeight">The total heights of the layers under this one. Basically this is the ground under the water.</param>
public override void DoUpdate(float dt, float dx, float[][] lowerLayersHeight) {
    _fluidTimer.Start();

    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 and Velocity Field Update
            // ----------------------------------------------------------------------------------------
            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);

    _fluidTimer.Stop();
}