def METHOD_NAME(particles):
    num = particles['rho'].shape[0]
    mass = dx*dy*particles['rho']
    position = particles["position"]
    velocity = particles["velocity"] 
    rij = position[:, np.newaxis, :] - position
    gk = gradkernel(rij, H)
    vij = velocity[:, np.newaxis, :] - velocity
    newrho = np.einsum('ijk, ijk, j, ij-> i', vij, rij, mass, gk) 
    particles["rho"] += dt * newrho