                             &rho, &u0, &u1, &u2);
 
                 // Compute equilibrium
-                double feq_ne, feq_e, feq_se, feq_n, feq_c, feq_s, feq_nw, feq_w, feq_sw, 
-                       feq_te, feq_tn, feq_tc, feq_ts, feq_tw, 
-                       feq_be, feq_bn, feq_bc, feq_bs, feq_bw;
                 d_equilibrium(&feq_ne, &feq_e, &feq_se, &feq_n, &feq_c, &feq_s, &feq_nw, &feq_w, &feq_sw, 
                               &feq_te, &feq_tn, &feq_tc, &feq_ts, &feq_tw, 
                               &feq_be, &feq_bn, &feq_bc, &feq_bs, &feq_bw, 
                               rho, u0, u1, u2);       
-
-                // Collision step
+                
+                d_vars.u.u0[gi] = u0;
+                d_vars.u.u1[gi] = u1;
+                d_vars.u.u2[gi] = u2;
+
+                // Collision
                 fout_ne = __dadd_rn(__dmul_rn(-OMEGA, __dadd_rn(fin_ne, - feq_ne)), fin_ne);
                 fout_e  = __dadd_rn(__dmul_rn(-OMEGA, __dadd_rn(fin_e , - feq_e )), fin_e );
                 fout_se = __dadd_rn(__dmul_rn(-OMEGA, __dadd_rn(fin_se, - feq_se)), fin_se);

    lbm_sim->dimComputationGrid = dimComputationGrid;
    lbm_sim->dimComputationBlock = dimComputationBlock;

    dim3 dimRightWallGrid(1, min((unsigned long)ny, (unsigned long)prop.maxGridSize[1]));
    dim3 dimRightWallBlock(1);
    lbm_sim->dimRightWallGrid = dimRightWallGrid;
    lbm_sim->dimRightWallBlock = dimRightWallBlock;

    lbm_sim->shared_mem_size = 6 * sizeof(double) * BLOCK_SIZE;


    return lbm_sim;
}

#define OMEGA    ((double)1. / (3*(NULB)+0.5))  // Relaxation parameter

#define BLOCK_SIZE 64

#ifndef SHARED_MEM_SIZE
#define SHARED_MEM_SIZE 0
#endif

#define SQUARE(a) ((a)*(a))
#define GPU_SQUARE(a) __dmul_rn(a,a)
