 #define NULB     ((ULB) * (R) / (RE))   // Viscoscity in lattice units
 #define OMEGA    ((double)1. / (3*(NULB)+0.5))  // Relaxation parameter
 
+#define NB_THREADS 100
+
 #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;


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

        // Compute macroscopic variables, density and velocity
        macroscopic(d_vars->fin[x][y], &d_vars->rho[x][y], d_vars->u[x][y]);


__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);

        // Compute equilibrium
        d_equilibrium(d_vars->feq[x][y], d_vars->rho[x][y], d_vars->u[x][y]);
    }
