I found the problem. My step was not working properly. new step function for Fluid:

fn step(&mut self) {
        const ITER: i32 = 16;

        diffuse(1, &mut self.vx0, &mut self.vx, self.visc, self.dt, ITER);
        diffuse(2, &mut self.vy0, &mut self.vy, self.visc, self.dt, ITER);

        project(
            &mut self.vx0,
            &mut self.vy0,
            &mut self.vx,
            &mut self.vy,
            ITER,
        );

        advect(
            1,
            &mut self.vx,
            Axis::X,
            &mut self.vx0,
            &mut self.vy0,
            self.dt,
        );
        advect(
            2,
            &mut self.vy,
            Axis::Y,
            &mut self.vx0,
            &mut self.vy0,
            self.dt,
        );

        project(
            &mut self.vx,
            &mut self.vy,
            &mut self.vx0,
            &mut self.vy0,
            ITER,
        );

        diffuse(0, &mut self.s, &mut self.density, self.diff, self.dt, ITER);
        advect2(
            0,
            &mut self.density,
            &mut self.s,
            &mut self.vx,
            &mut self.vy,
            self.dt,
        );

        set_bnd(1, &mut self.vx);
        set_bnd(2, &mut self.vy);
        set_bnd(0, &mut self.density);
    }

and added an advect2 function which functions same way as defined in the jos stam solver on which mine is based on. Here is the code:

fn advect2<'a>(
    b: usize,
    d: &mut Array2D,
    d0: &mut Array2D,
    vx: &'a mut Array2D,
    vy: &'a mut Array2D,
    dt: f32,
) {
    let dtx = dt * (N - 2) as f32;
    let dty = dt * (N - 2) as f32;

    let n_float = N as f32;

    let (mut i0, mut i1, mut j0, mut j1);
    let (mut tmp1, mut tmp2, mut x, mut y);
    let (mut s0, mut s1, mut t0, mut t1);

    for i in 1..(N - 1) {
        for j in 1..(N - 1) {
            tmp1 = dtx * vx[i][j];
            tmp2 = dty * vy[i][j];

            x = i as f32 - tmp1;
            y = j as f32 - tmp2;

            x = clamp(x, 0.5, n_float + 0.5);
            i0 = x.floor();
            i1 = i0 + 1.0;

            y = clamp(y, 0.5, n_float + 0.5);
            j0 = y.floor();
            j1 = j0 + 1.0;

            s1 = x - i0;
            s0 = 1.0 - s1;
            t1 = y - j0;
            t0 = 1.0 - t1;

            let i0i = i0 as usize;
            let i1i = i1 as usize;
            let j0i = j0 as usize;
            let j1i = j1 as usize;

            d[i][j] = s0 * (t0 * d0[i0i][j0i] + t1 * d0[i0i][j1i])
                + s1 * (t0 * d0[i1i][j0i] + t1 * d0[i1i][j1i]);
        }
    }

    set_bnd(b, d);
}