code: def advect(f, vx, vy):
      
      rows, cols = f.shape
      cell_xs, cell_ys = np.meshgrid(np.arange(cols), np.arange(rows))
      center_xs = (cell_xs - vx).ravel()
      center_ys = (cell_ys - vy).ravel()
  
      # Compute indices of source cells.
      left_ix = np.floor(center_ys).astype(int)
      top_ix  = np.floor(center_xs).astype(int)
      rw = center_ys - left_ix              # Relative weight of right-hand cells.
      bw = center_xs - top_ix               # Relative weight of bottom cells.
      left_ix  = np.mod(left_ix,     rows)  # Wrap around edges of simulation.
      right_ix = np.mod(left_ix + 1, rows)
      top_ix   = np.mod(top_ix,      cols)
      bot_ix   = np.mod(top_ix  + 1, cols)
  
      # A linearly-weighted sum of the 4 surrounding cells.
      flat_f = (1 - rw) * ((1 - bw)*f[left_ix,  top_ix] + bw*f[left_ix,  bot_ix]) \
                   + rw * ((1 - bw)*f[right_ix, top_ix] + bw*f[right_ix, bot_ix])
      return np.reshape(flat_f, (rows, cols))
docs: Move field f according to x and y velocities (u and v)
         using an implicit Euler integrator.