     HANDLE_ERROR(cudaMemcpy(d_vars, h_vars, sizeof(lbm_vars), cudaMemcpyHostToDevice));
 
     pgm_image* pgm = pgm_create(NX, NY);
+
+#ifndef COMPUTE_ON_CPU
+    dim3 dimBlock(1);
+    dim3 dimGrid(getThreads(NX, NY));
+#endif
 
     for (int time = 0; time < max_iter; time++) {
         RUN_KERNEL_1D(lbm_right_wall,             NY, d_vars);

    u[1] /= *rho;
}

#ifndef COMPUTE_ON_CPU
__global__ void lbm_right_wall(lbm_vars *d_vars)
#else
void lbm_right_wall(lbm_vars *d_vars, int y)
#endif
{
#ifndef COMPUTE_ON_CPU
    int y = threadIdx.x;
#endif

    // Right wall: outflow condition.
    for (int i = 0; i < 3; i++) {
        int f = d_consts.col[2][i];
        d_vars->fin[NX-1][y][f] = d_vars->fin[NX-2][y][f]; // TODO Test: retirer condition bord
    }
}

#ifndef COMPUTE_ON_CPU
__global__ void lbm_macro_and_left_wall(lbm_vars *d_vars)
#else
void lbm_macro_and_left_wall(lbm_vars *d_vars, int x, int y)
#endif
{
#ifndef COMPUTE_ON_CPU
    int y = threadIdx.x / NX;
    int x = threadIdx.x % NX;
#endif

    // Compute macroscopic variables, density and velocity
    macroscopic(d_vars->fin[x][y], &d_vars->rho[x][y], d_vars->u[x][y]);
    
    // Left wall: inflow condition
    for (size_t d = 0; d < 2; d++) {
        d_vars->u[0][y][d] = d_vars->vel[0][y][d]; // TODO: collision!
    }
 }

#ifndef COMPUTE_ON_CPU
__global__ void lbm_density(lbm_vars *d_vars)
#else
void lbm_density(lbm_vars *d_vars, int y)
#endif
{
#ifndef COMPUTE_ON_CPU
    int y = threadIdx.x;
#endif
    // Calculate the density
    double s2 = 0, s3 = 0;

        d_vars->fin[x_dst][y_dst][f] = d_vars->fout[x][y][f];
    }

#ifndef COMPUTE_ON_CPU
    }
#endif
}

void print_variables(lbm_vars *d_vars, lbm_vars *h_vars, double var[NX][NY][9]) {
