ABM and CellListMap

https://github.com/m3g/CellListMap.jl can quickly calculate partical pairs interactions within a cutoff.

Modeling the particle repulsion:

\begin{align} U(r) &= k_i k_j (r^2 - (r_i - r_j)^2)^2, \quad \text{for} r ≤ (r_i + r_j) \\ U(r) &= 0, \quad \text{otherwise} \end{align}

(1)

using Agents
using CellListMap
using Base64
using CairoMakie
CairoMakie.activate!(px_per_unit = 1.0)

The helper function displays video files in Jupyter notebooks

function display_mp4(filename)
    display("text/html", string("""<video autoplay controls><source src="data:video/x-m4v;base64,""",
        base64encode(open(read, filename)),
        """" type="video/mp4"></video>"""))
end

display_mp4 (generic function with 1 method)

Define particle agents

@agent struct Particle(ContinuousAgent{2,Float64})
    r::Float64 ## radius
    k::Float64 ## repulsion force constant
    mass::Float64
end

Particle(; vel, r, k, mass) = (vel, r, k, mass)

Main.var"##234".Particle

Building the model

function initialize_bouncing(;
    number_of_particles=10_000,
    sides=SVector(500.0, 500.0),
    dt=0.001,
    max_radius=10.0,
    parallel=true
)
    # initial random positions
    positions = [sides .* rand(SVector{2,Float64}) for _ in 1:number_of_particles]

    # We will use CellListMap to compute forces, with similar structure as the positions
    forces = similar(positions)

    # Space for the agents
    space2d = ContinuousSpace(sides; periodic=true)

    # Initialize CellListMap particle system
    system = ParticleSystem(
        positions=positions,
        unitcell=sides,
        cutoff=2 * max_radius,
        output=forces,
        output_name=:forces, ## allows the system.forces alias for clarity
        parallel=parallel,
    )

    # define the ABModel properties
    # The system field contains the data required for CellListMap.jl
    properties = (;dt, number_of_particles, system)

    model = StandardABM(Particle,
        space2d;
        agent_step!,
        model_step!,
        agents_first = false,
        properties=properties
    )

    # Create active agents
    for id in 1:number_of_particles
        pos = positions[id]
        prop_particle = Particle(
            r = (0.5 + 0.9 * rand()) * max_radius,
            k = 10 + 20 * rand(), ## random force constants
            mass = 10.0 + 100 * rand(), ## random masses
            vel = 100 * randn(SVector{2}) ## initial velocities)
            )
        add_agent!(pos, Particle, model, prop_particle...)
    end

    return model
end

initialize_bouncing (generic function with 1 method)

Computing the repulsion force It must follow CellListMap API and return a array for forces acting upon the particles

function calc_forces!(x, y, i, j, d2, forces, model)
    p_i = model[i]
    p_j = model[j]
    d = sqrt(d2)
    if d ≤ (p_i.r + p_j.r)
        dr = y - x ## x and y are minimum-image relative coordinates
        fij = 2 * (p_i.k * p_j.k) * (d2 - (p_i.r + p_j.r)^2) * (dr / d)
        forces[i] += fij
        forces[j] -= fij
    end
    return forces
end

calc_forces! (generic function with 1 method)

Update the pairwise forces using CellListMap API

function model_step!(model::ABM)
    map_pairwise!(
        (x, y, i, j, d2, forces) -> calc_forces!(x, y, i, j, d2, forces, model),
        model.system,
    )
    return nothing
end

model_step! (generic function with 1 method)

Update agent positions and velocity

function agent_step!(agent, model::ABM)
    id = agent.id
    dt = abmproperties(model).dt
    force = model.system.forces[id]
    acc = force / agent.mass
    # Update positions and velocities
    vel = agent.vel + acc * dt
    x = agent.pos + vel * dt + (acc / 2) * dt^2
    x = normalize_position(x, model)  ## Wraps agent position
    agent.vel = vel
    move_agent!(agent, x, model)
    # !!! Remember to update positions in the ParticleSystem
    model.system.positions[id] = agent.pos
    return nothing
end

agent_step! (generic function with 1 method)

Run simulation

function simulate(model=nothing; nsteps=1_000, number_of_particles=10_000)
    if isnothing(model)
        model = initialize_bouncing(number_of_particles=number_of_particles)
    end
    Agents.step!(model, nsteps)
end

simulate (generic function with 2 methods)

Test the performance

model = initialize_bouncing(number_of_particles=10_000)
@time simulate(model)

 10.880838 seconds (317.13 k allocations: 26.921 MiB, 3.09% compilation time)

StandardABM with 10000 agents of type Particle
 agents container: Dict
 space: periodic continuous space with [500.0, 500.0] extent and spacing=25.0
 scheduler: fastest
 properties: dt, number_of_particles, system

Visualize

model = initialize_bouncing(number_of_particles=1000)
abmvideo(
    "cellistmap.mp4",
    model;
    framerate=20, frames=200, dt=5,
    title="Softly bouncing particles",
    agent_size=p -> p.r,
    agent_color=p -> p.k
)

display_mp4("cellistmap.mp4")

This notebook was generated using Literate.jl.