                 d_vars->fout[x][y][f] = __dadd_rn(__dmul_rn(-OMEGA, __dadd_rn(d_vars->fin[x][y][f], - feq[f])), d_vars->fin[x][y][f]);
             }
         }
+
+		d_vars->u[x][y][0] = u0;
+		d_vars->u[x][y][1] = u1;
     }
 }
 

        }

        // Compute macroscopic variables, density and velocity
        macroscopic(d_vars->fin[x][y], &d_vars->rho[x][y], d_vars->u[x][y]);
        
        if (x == 0) {
            // Left wall: inflow condition
            for (size_t d = 0; d < 2; d++) {
                d_vars->u[0][y][d] = d_vars->vel[0][y][d];
            }   

            // Calculate the density
            double s2 = 0, s3 = 0;

                s2 += d_vars->fin[0][y][d_consts.col[1][i]];
                s3 += d_vars->fin[0][y][d_consts.col[2][i]];
            }
            d_vars->rho[0][y] = 1./(1 - u0) * (s2 + 2*s3);
        }

        // Compute equilibrium
        double feq[9];
        d_equilibrium(feq, d_vars->rho[x][y], u0, u1);

        if (x == 0) {
            for (size_t i = 0, f = d_consts.col[0][i]; i < 3; f = d_consts.col[0][++i]) {
