Verlet list

.gitignore

@@ -0,0 +1,2 @@
+Manifest.toml
+exports

Project.toml

@@ -0,0 +1,11 @@
+[deps]
+BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
+Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
+GLMakie = "e9467ef8-e4e7-5192-8a1a-b1aee30e663a"
+JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
+LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
+ProfileView = "c46f51b8-102a-5cf2-8d2c-8597cb0e0da7"
+ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca"
+Revise = "295af30f-e4ad-537b-8983-00126c2a3abe"
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
+Traceur = "37b6cedf-1f77-55f8-9503-c64b63398394"

src/Particle.jl

@@ -0,0 +1,64 @@
+mutable struct Particle
+    id::Int64
+    c::Vector{Float64} # Center
+    tmp_c::Vector{Float64} # Temporary center
+
+    φ::Float64 # Angle
+
+    function Particle(; id::Int64, c::SVector{2,Float64}, φ::Float64)
+        return new(id, c, c, φ)
+    end
+end
+
+get_c(p::Particle) = p.c
+get_c_x(p::Particle) = p.c[1]
+get_c_y(p::Particle) = p.c[2]
+get_φ(p::Particle) = p.φ
+
+function restrict_coordinate(value::Float64; l::Float64)
+    if value < -l
+        value += 2 * l
+    elseif value >= l
+        value -= 2 * l
+    end
+
+    return value
+end
+
+function restrict_coordinates!(p::Particle; l::Float64)
+    @simd for i in 1:2
+        p.tmp_c[i] = restrict_coordinate(p.tmp_c[i]; l=l)
+    end
+
+    return nothing
+end
+
+function minimum_image_coordinate(value::Float64; l::Float64)
+    if value <= -l
+        value += 2 * l
+    elseif value > l
+        value -= 2 * l
+    end
+
+    return value
+end
+
+function minimum_image(v::SVector{2,Float64}; l::Float64)
+    return minimum_image_coordinate.(v; l=l)
+end
+
+function are_overlapping(p1::Particle, p2::Particle, overlapping_r²::Float64, l::Float64)
+    r⃗₁₂ = SVector{2}(p2.c) - SVector{2}(p1.c) # 1 -> 2
+
+    r⃗₁₂ = minimum_image(r⃗₁₂; l=l)
+
+    distance² = r⃗₁₂[1]^2 + r⃗₁₂[2]^2
+
+    return (distance² < overlapping_r², r⃗₁₂, distance²)
+end
+
+function update!(p::Particle)
+    @simd for i in 1:2
+        p.c[i] = p.tmp_c[i]
+    end
+end

src/PreVector.jl

@@ -0,0 +1,19 @@
+mutable struct PreVector{T}
+    last_ind::UInt64
+    v::Vector{T}
+
+    PreVector(T, N) = new{T}(UInt64(0), Vector{T}(undef, N))
+end
+
+function push!(pv::PreVector{T}, x::T) where {T}
+    pv.last_ind += 1
+    pv.v[pv.last_ind] = x
+
+    return nothing
+end
+
+function reset!(pv::PreVector{T}) where {T}
+    pv.last_ind = 0
+
+    return nothing
+end

src/ReCo.jl

@@ -0,0 +1,24 @@
+using Pkg
+Pkg.activate(".")
+
+using Random, Distributions
+using JLD2
+using ProgressMeter
+using GLMakie
+using LaTeXStrings
+using StaticArrays
+
+import Dates: now, CompoundPeriod, canonicalize
+import Base.push!
+
+# Development deps
+using Revise
+using BenchmarkTools
+
+includet("PreVector.jl")
+includet("Particle.jl")
+includet("setup.jl")
+includet("simulation.jl")
+includet("data.jl")
+includet("animation.jl")
+includet("run.jl")

src/animation.jl

@@ -0,0 +1,72 @@
+function animate(sol::Solution, args, name_part::String; framerate::Int64=10)
+    println("Generating animation...")
+
+    fig = Figure()
+    ax = Axis(
+        fig[1, 1];
+        limits=(-args.l, args.l, -args.l, args.l),
+        aspect=AxisAspect(1),
+        xlabel=L"x",
+        ylabel=L"y",
+    )
+
+    animation_path = "exports/$name_part.mkv"
+
+    record(fig, animation_path; framerate=framerate) do io
+        old_cx = zeros(args.N)
+        old_cy = zeros(args.N)
+
+        new_cx = zeros(args.N)
+        new_cy = zeros(args.N)
+
+        segments_x = zeros(2 * args.N)
+        segments_y = zeros(2 * args.N)
+
+        done_first_iter = false
+
+        @showprogress 0.5 for frame in 1:args.n_frames
+            new_cx .= get_c_x.(sol.position[frame])
+            new_cy .= get_c_y.(sol.position[frame])
+
+            if done_first_iter
+                empty!(ax)
+
+                for i in 1:(args.N)
+                    segments_x[2 * i - 1] = old_cx[i]
+                    segments_x[2 * i] = new_cx[i]
+
+                    segments_y[2 * i - 1] = old_cy[i]
+                    segments_y[2 * i] = new_cy[i]
+                end
+
+                linesegments!(ax, segments_x, segments_y; color=1:(args.N))
+            else
+                done_first_iter = true
+                println("Started recording!")
+            end
+
+            old_cx .= new_cx
+            old_cy .= new_cy
+
+            poly!(
+                ax,
+                [Circle(Point2(x, y), args.particle_diameter / 2) for (x, y) in zip(new_cx, new_cy)];
+                color=1:(args.N),
+            )
+
+            # DEBUG BEGIN
+            for (x, y) in zip(new_cx, new_cy)
+                arc!(ax, Point2(x, y), args.interaction_r, 0.0, 2 * π)
+            end
+            # DEBUG END
+
+            ax.title = "t = $(round(sol.t[frame], digits=2))"
+
+            recordframe!(io)
+        end
+    end
+
+    println("Animation done and saved to $animation_path.")
+
+    return nothing
+end

src/benchmark.jl

@@ -0,0 +1 @@
+using BenchmarkTools

src/data.jl

@@ -0,0 +1,20 @@
+struct Solution
+    t::Vector{Float64}
+    position::Matrix{Vector{Float64}}
+    φ::Matrix{Float64}
+
+    function Solution(N::Int64, n_frames::Int64)
+        return new(zeros(n_frames), [zeros(2) for i in 1:N, j in 1:n_frames], zeros((N, n_frames)))
+    end
+end
+
+function save_data_jld(sol, args, name_part)
+    println("Saving data...")
+
+    data_path = "exports/$name_part.jld2"
+    jldsave(data_path; sol, args, name_part)
+
+    println("Data saved to $data_path.")
+
+    return nothing
+end

src/run.jl

@@ -0,0 +1,63 @@
+function run(;
+    N::Int64,
+    T::Float64,
+    v::Float64=20.0,
+    δt::Float64=2e-5,
+    save_at::Float64=0.1,
+    framerate::Int64=0,
+    save_data::Bool=false,
+    n_steps_before_verlet_list_update::Int64=100,
+    )
+
+    Random.seed!(42)
+
+    μ = 1.0
+    D₀ = 1.0
+    particle_diameter = 1.0
+    Dᵣ = 3 * D₀ / (particle_diameter^2)
+    σ = 1.0
+    ϵ = 100.0
+    interaction_r = 2^(1 / 6) * σ
+
+    grid_n = round(Int64, ceil(sqrt(N)))
+
+    N = grid_n^2
+
+    l = float(sqrt(N))
+    grid_box_width = 2 * l / grid_n
+
+    skin_r = 1.1 * (2 * v * n_steps_before_verlet_list_update * δt + 2 * interaction_r)
+
+    args = (
+        v = v,
+        c₁ = 4 * ϵ * 6 * σ^6 * δt * μ,
+        c₂ = 2 * σ^6,
+        c₃ = sqrt(2 * D₀ * δt),
+        c₄ = sqrt(2 * Dᵣ * δt),
+        vδt = v * δt,
+        μ = μ,
+        interaction_r = interaction_r,
+        interaction_r² = interaction_r^2,
+        N = N,
+        l = l,
+        particle_diameter = particle_diameter,
+        particles = generate_particles(grid_n, grid_box_width, l),
+        skin_r² = skin_r^2,
+        verlet_list = [PreVector(Int64, N - 1) for i in 1:N],
+        n_frames = floor(Int64, T / save_at) + 1,
+    )
+
+    sol, end_time = simulate(args, save_at, δt, T, n_steps_before_verlet_list_update)
+
+    name_part = "$(end_time)_T=$(T)_N=$(N)_v=$(v)_dt=$(δt)"
+
+    if save_data
+        save_data_jld(sol, args, name_part)
+    end
+
+    if framerate > 0
+        generate_animation(sol, args, name_part; framerate=framerate)
+    end
+
+    return (sol = sol, args = args, name_part = name_part)
+end

src/setup.jl

@@ -0,0 +1,27 @@
+function initial_particle_grid_pos(i, j; grid_box_width, l)
+    dim1_pos(x) = (x - 0.5) * grid_box_width - l
+
+    return dim1_pos.(SVector(i, j))
+end
+
+function generate_particles(grid_n, grid_box_width, l)
+    particles = Vector{Particle}(undef, grid_n^2)
+
+    particle_pos_in_grid_dim(i) = (i - 0.5) * grid_box_width - l
+
+    id = 1
+
+    for i in 1:grid_n
+        for j in 1:grid_n
+            particles[id] = Particle(;
+                id=id,
+                c=initial_particle_grid_pos(i, j; grid_box_width=grid_box_width, l=l),
+                φ=rand(Uniform(-π, π)),
+            )
+
+            id += 1
+        end
+    end
+
+    return particles
+end

src/simulation.jl

@@ -0,0 +1,108 @@
+rand_normal01() = rand(Normal(0, 1))
+
+function update_verlet_list!(args)
+    @simd for pv in args.verlet_list
+        reset!(pv)
+    end
+
+    for i in 1:(args.N - 1)
+        for j in (i + 1):args.N
+            p1 = args.particles[i]
+            p2 = args.particles[j]
+
+            overlapping, r⃗₁₂, distance² = are_overlapping(p1, p2, args.skin_r², args.l)
+
+            if overlapping
+                push!(args.verlet_list[i], j)
+                push!(args.verlet_list[j], i)
+            end
+        end
+    end
+end
+
+function euler!(args)
+    Threads.@threads for p in args.particles
+        verlet_list = args.verlet_list[p.id]
+
+        for i in 1:verlet_list.last_ind
+            p2 = args.particles[verlet_list.v[i]]
+
+            overlapping, r⃗₁₂, distance² = are_overlapping(p, p2, args.interaction_r², args.l)
+
+            if overlapping
+                c = args.c₁ / (distance²^4) * (args.c₂ / (distance²^3) - 1)
+                @simd for j in 1:2
+                    p.tmp_c[j] -= c * r⃗₁₂[j]
+                end
+            end
+        end
+
+        e = SVector(cos(p.φ), sin(p.φ))
+
+        @simd for i in 1:2
+            p.tmp_c[i] += args.vδt * e[i] + args.c₃ * rand_normal01()
+        end
+
+        p.φ += args.c₄ * rand_normal01()
+
+        restrict_coordinates!(p; l=args.l)
+    end
+
+    @simd for particle in args.particles
+        update!(particle)
+    end
+
+    return nothing
+end
+
+function simulate(args, save_at::Float64, δt::Float64, T::Float64, n_steps_before_verlet_list_update::Int64)
+    sol = Solution(args.N, args.n_frames)
+    save_timer = save_at
+    frame = 0
+
+    start_time = now()
+    println("Started simulation at $start_time.")
+
+    update_verlet_list_at = n_steps_before_verlet_list_update * δt
+    update_verlet_list_timer = update_verlet_list_at
+
+    @showprogress 0.2 for t in 0:δt:T
+        if save_timer >= save_at
+            frame += 1
+
+            sol.t[frame] = t
+
+            @simd for p in args.particles
+                p_id = p.id
+
+                pre_pos = sol.position[p_id, frame]
+                
+                @simd for i in 1:2
+                    pre_pos[i] = p.c[i]
+                end
+
+                sol.φ[p_id, frame] = p.φ
+            end
+
+            save_timer = 0.0
+        else
+            save_timer += δt
+        end
+
+        if update_verlet_list_timer >= update_verlet_list_at
+            update_verlet_list!(args)
+            
+            update_verlet_list_timer = 0.0
+        else
+            update_verlet_list_timer += δt
+        end
+
+        euler!(args)
+    end
+
+    end_time = now()
+    elapsed_time = canonicalize(CompoundPeriod(end_time - start_time))
+    println("Done simulation at $end_time and took $elapsed_time.")
+
+    return (sol, end_time)
+end