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

    *u2 = (te + tn + tc + ts + tw - be - bn - bc - bs - bw) / *rho;
}

__global__ void lbm_right_wall(lbm_lattices* f)
{
    for (int z = blockIdx.z; z < NZ; z+=gridDim.z) {
        for (int y = blockIdx.y; y < NY; y+=gridDim.y) {
           // Right wall: outflow condition.
           f->nw[IDX(NX-1,y,z)] = f->nw[IDX(NX-2,y,z)];
           f->w [IDX(NX-1,y,z)] = f->w [IDX(NX-2,y,z)];
           f->sw[IDX(NX-1,y,z)] = f->sw[IDX(NX-2,y,z)];
        }
    }
}

__global__ void lbm_computation(lbm_vars *d_vars, lbm_lattices* f0, lbm_lattices* f1)
{
    int tix = threadIdx.x;

                            fin_be, fin_bn, fin_bc, fin_bs, fin_bw,
                            &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, 
