 void lbm_simulation_update(lbm_simulation* lbm_sim)
 {
     if (lbm_sim->switch_f0_f1) {
-//        HANDLE_KERNEL_ERROR(lbm_right_wall<<<lbm_sim->dimRightWallGrid, lbm_sim->dimRightWallBlock>>>(&lbm_sim->d_vars->f1));
         HANDLE_KERNEL_ERROR(lbm_computation<<<lbm_sim->dimComputationGrid, lbm_sim->dimComputationBlock, lbm_sim->shared_mem_size>>>(lbm_sim->d_vars, &lbm_sim->d_vars->f1, &lbm_sim->d_vars->f0));
     } else {
-//        HANDLE_KERNEL_ERROR(lbm_right_wall<<<lbm_sim->dimRightWallGrid, lbm_sim->dimRightWallBlock>>>(&lbm_sim->d_vars->f0));
         HANDLE_KERNEL_ERROR(lbm_computation<<<lbm_sim->dimComputationGrid, lbm_sim->dimComputationBlock, lbm_sim->shared_mem_size>>>(lbm_sim->d_vars, &lbm_sim->d_vars->f0, &lbm_sim->d_vars->f1));
     }
 

                            fin_be, fin_bn, fin_bc, fin_bs, fin_bw,
                            &rho, &u0, &u1, &u2);

                
//                if (x == 0) {
//                    // Left wall: inflow condition
//                    u0 = d_consts.vel[y];
//                    u1 = 0;
//                    u2 = 0;
//
//                    // Calculate the density
//                    double s2 = fin_n  + fin_c + fin_s + fin_tn + fin_tc + fin_ts + fin_bn + fin_bc + fin_bs;
//                    double s3 = fin_nw + fin_w + fin_sw + fin_tw + fin_bw;
//                    rho = 1./(1 - u0) * (s2 + 2*s3);
//                }
//          

                // 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_te, &feq_tn, &feq_tc, &feq_ts, &feq_tw, 
                              &feq_be, &feq_bn, &feq_bc, &feq_bs, &feq_bw, 
                              rho, u0, u1, u2);       

                if (d_vars->obstacles[IDX(x, y, z)]) {
                    // Bounce-back condition for obstacle
