docs: def project(vx, vy):
      """Project the velocity field to be approximately mass-conserving,
         using a few iterations of Gauss-Seidel."""
      p = np.zeros(vx.shape)
      h = 1.0/vx.shape[0]
      div = -0.5 * h * (np.roll(vx, -1, axis=0) - np.roll(vx, 1, axis=0)
                      + np.roll(vy, -1, axis=1) - np.roll(vy, 1, axis=1))
  
      for k in range(10):
          p = (div + np.roll(p, 1, axis=0) + np.roll(p, -1, axis=0)
                   + np.roll(p, 1, axis=1) + np.roll(p, -1, axis=1))/4.0
  
      vx -= 0.5*(np.roll(p, -1, axis=0) - np.roll(p, 1, axis=0))/h
      vy -= 0.5*(np.roll(p, -1, axis=1) - np.roll(p, 1, axis=1))/h
      return vx, vy