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

struct lbm_simulation{
    lbm_vars h_vars, d_vars;
    dim3 dimComputationGrid, dimComputationBlock;
    dim3 dimRightWallGrid, dimRightWallBlock;
    size_t shared_mem_size;
    bool switch_f0_f1;
    size_t nx, ny, nz;

                f1.bc[IDX(x  , y  , z-1, nx,ny,nz)] = fout_bc;
                f1.bs[IDX(x  , y-1, z-1, nx,ny,nz)] = fout_bs;
                f1.bw[IDX(x-1, y  , z-1, nx,ny,nz)] = fout_bw;
           }
        }
    }
