 #define SQUARE(a) ((a)*(a))
 #define GPU_SQUARE(a) (__dmul_rn(a,a))
 
+#define NB_THREADS 100
+
+#define INDEX_2D_FROM_1D(x, y, i) (y) = (i)/(NX), (x) = (i)%(NX)
+
 typedef enum { OUT_FIN, OUT_IMG, OUT_UNK } out_mode;
 
 typedef struct {

    lbm_vars *d_vars;
    HANDLE_ERROR(cudaMalloc(&d_vars, sizeof(lbm_vars)));
    HANDLE_ERROR(cudaMemcpy(d_vars, h_vars, sizeof(lbm_vars), cudaMemcpyHostToDevice));
    for (int time = 0; time < max_iter; time++) {
        HANDLE_KERNEL_ERROR(         lbm_right_wall<<<1,   NY>>>(d_vars));        
        HANDLE_KERNEL_ERROR(lbm_macro_and_left_wall<<<1,NX*NY>>>(d_vars));
        HANDLE_KERNEL_ERROR(            lbm_density<<<1,   NY>>>(d_vars));
        HANDLE_KERNEL_ERROR(      lbm_equilibrium_1<<<1,NX*NY>>>(d_vars));
        HANDLE_KERNEL_ERROR(      lbm_equilibrium_2<<<1,   NY>>>(d_vars));
        HANDLE_KERNEL_ERROR(          lbm_collision<<<1,NX*NY>>>(d_vars));
        HANDLE_KERNEL_ERROR(          lbm_streaming<<<1,NX*NY>>>(d_vars));
    }

    if (out == OUT_FIN) {

}

__device__ static void macroscopic(double* fin, double* rho, double* u)
{   
    *rho = u[0] = u[1] = 0;

    for (int f = 0; f < 9; f++) {
