 
 __global__ void lbm_right_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);
-        
+    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) {
             // Right wall: outflow condition.
             for (int i = 0; i < 3; i++) {


__global__ void lbm_equilibrium_2(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);
   
        if (x == 0) {
            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]];

#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 INDEX_2D_FROM_1D(x, y, i) (y) = (i)/(NX), (x) = (i)%(NX)

typedef enum { OUT_FIN, OUT_IMG, OUT_UNK } out_mode;
