                 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 - d_vars->u[0][y][0]) * (s2 + 2*s3);
+            d_vars->rho[0][y] = 1./(1 - u0) * (s2 + 2*s3);
         }
 
         // Compute equilibrium
         double feq[9];
-        d_equilibrium(feq, d_vars->rho[x][y], d_vars->u[x][y]);
+        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]) {

    }                                                                
}

__device__ static void d_equilibrium(double* feq, double rho, double* u)
{
    double usqr = __dmul_rn(3./2, __dadd_rn( GPU_SQUARE(u[0]), GPU_SQUARE(u[1]) ));
                                                                     
    for (int f = 0; f < 9; f++) {
        double cu = 3 * ( d_consts.v[0][f] * u[0] + d_consts.v[1][f] * u[1] );
        feq[f] = rho * d_consts.t[f] * ( 1 + cu + 0.5 * SQUARE(cu) - usqr );
    }                                                                
}

__device__ static void macroscopic(double* fin, double* rho, double* u)
{   
    *rho = u[0] = u[1] = 0;

    for (int f = 0; f < 9; f++) {
        *rho += fin[f];

        u[0] += d_consts.v[0][f] * fin[f];
        u[1] += d_consts.v[1][f] * fin[f];
    }
    
    u[0] /= *rho;
    u[1] /= *rho;
}

__global__ void lbm_computation(lbm_vars *d_vars)

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

