 __global__ void lbm_streaming(lbm_vars *d_vars)
 {
 
-    for (int i = threadIdx.x + blockIdx.x * blockDim.x; i < NX*NY; i += blockDim.x * gridDim.x) {
-        int x, y;
-        INDEX_2D_FROM_1D(x, y, i);
+    for (int idx = threadIdx.x + blockIdx.x * blockDim.x; idx < NX*NY; idx += blockDim.x * gridDim.x) {
+        int x, y;
+        INDEX_2D_FROM_1D(x, y, idx);
 
         // Streaming step
         for (size_t f = 0; f < 9; f++) {


__global__ void lbm_equilibrium_1(lbm_vars *d_vars)
{
    for (int i = threadIdx.x + blockIdx.x * blockDim.x; i < NX*NY; i += blockDim.x * gridDim.x) {
        int x, y;
        INDEX_2D_FROM_1D(x, y, i);
   
        // Compute equilibrium
        d_equilibrium(d_vars->feq[x][y], d_vars->rho[x][y], d_vars->u[x][y]);
    }


__global__ void lbm_density(lbm_vars *d_vars)
{
    for (int idx = threadIdx.x + blockIdx.x * blockDim.x; idx < NX*NY; idx += blockDim.x * gridDim.x) {
        int x, y;
        INDEX_2D_FROM_1D(x, y, idx);

        if (x == 0) {
            // Calculate the density
            double s2 = 0, s3 = 0;
