     }
 }
 
+int getThreads(int width, int height) {
+    int dev, threads;
+    cudaDeviceProp prop;
+    HANDLE_ERROR( cudaGetDevice(&dev) );
+    HANDLE_ERROR( cudaGetDeviceProperties(&prop, dev) );
+
+    int maxThreads = min (prop.maxThreadsDim[0], prop.maxThreadsPerBlock);
+#ifdef NB_THREADS
+    threads = NB_THREADS;
+#else
+    threads = prop.maxThreadsDim[0];
+#endif
+
+    if (threads > maxThreads)
+        threads = maxThreads;
+
+    return min(threads, width*height);
+}
 
 int main(int argc, char * const argv[])
 {

    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++) {
