 
 __global__ void lbm_collision(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);
+
         for (size_t f = 0; f < 9; f++) {
             if (d_vars->obstacles[x][y]) {
                 // Bounce-back condition for obstacle


__global__ void lbm_equilibrium_2(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);
   
        if (x == 0) {
            for (size_t i = 0, f = d_consts.col[0][i]; i < 3; f = d_consts.col[0][++i]) {
                d_vars->fin[0][y][f] = d_vars->feq[0][y][f] + d_vars->fin[0][y][d_consts.opp[f]] - d_vars->feq[0][y][d_consts.opp[f]];

__global__ void lbm_streaming(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);

        // Streaming step
        for (size_t f = 0; f < 9; f++) {
