     HANDLE_ERROR( cudaDeviceSynchronize() ); \
 } while(0)
 
-#define RUN_KERNEL_1D(kernel, th1, ...) HANDLE_KERNEL_ERROR( kernel<<<1, th1>>>(__VA_ARGS__) )
-#define RUN_KERNEL_2D(kernel, th1, th2, ...) HANDLE_KERNEL_ERROR( kernel<<<1, th1*th2>>>(__VA_ARGS__) )
+#define RUN_KERNEL_1D(kernel, th1, ...) HANDLE_KERNEL_ERROR( kernel<<<dimBlock, dimGrid>>>(__VA_ARGS__) )
+#define RUN_KERNEL_2D(kernel, th1, th2, ...) HANDLE_KERNEL_ERROR( kernel<<<dimBlock, dimGrid>>>(__VA_ARGS__) )
 #endif
 
 // Constants

    HANDLE_ERROR(cudaMemcpy(d_vars, h_vars, sizeof(lbm_vars), cudaMemcpyHostToDevice));

    pgm_image* pgm = pgm_create(NX, NY);

    for (int time = 0; time < max_iter; time++) {
        RUN_KERNEL_1D(lbm_right_wall,             NY, d_vars);

        s3 += d_vars->fin[0][y][d_consts.col[2][i]];
    }
    d_vars->rho[0][y] = 1./(1 - d_vars->u[0][y][0]) * (s2 + 2*s3);

#ifndef COMPUTE_ON_CPU
        }   
    }
#endif
}

#ifdef COMPUTE_ON_CPU
void lbm_equilibrium_1(lbm_vars *d_vars, int x, int y)
{
#else
__global__ void lbm_equilibrium_1(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);
#endif
   
    // Compute equilibrium
    d_equilibrium(d_vars->feq[x][y], d_vars->rho[x][y], d_vars->u[x][y]);

#ifndef COMPUTE_ON_CPU
    }
#endif
}

#ifdef COMPUTE_ON_CPU
void lbm_equilibrium_2(lbm_vars *d_vars, int y)
{
#else 
__global__ void lbm_equilibrium_2(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) {
#endif
    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]];
    }
#ifndef COMPUTE_ON_CPU
        }
    }
#endif
}

#ifdef COMPUTE_ON_CPU
void lbm_collision(lbm_vars *d_vars, int x, int y)
{
#else
__global__ void lbm_collision(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);
#endif
    
    for (size_t f = 0; f < 9; f++) {
