diff --git a/Project.toml b/Project.toml index fe63a2f..96d8ebc 100644 --- a/Project.toml +++ b/Project.toml @@ -22,7 +22,6 @@ ProfileView = "c46f51b8-102a-5cf2-8d2c-8597cb0e0da7" ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" ReinforcementLearning = "158674fc-8238-5cab-b5ba-03dfc80d1318" -Revise = "295af30f-e4ad-537b-8983-00126c2a3abe" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" diff --git a/src/Particle.jl b/src/Particle.jl index 712b420..4a3f6da 100644 --- a/src/Particle.jl +++ b/src/Particle.jl @@ -13,46 +13,54 @@ mutable struct Particle end end -function restrict_coordinate(value::Float64, l::Float64) - if value < -l - value += 2 * l - elseif value >= l - value -= 2 * l +function restrict_coordinate(value::Float64, half_box_len::Float64) + if value < -half_box_len + value += 2 * half_box_len + elseif value >= half_box_len + value -= 2 * half_box_len end return value end -function restrict_coordinates!(v::SVector{2,Float64}, l::Float64) - return SVector(restrict_coordinate(v[1], l), restrict_coordinate(v[2], l)) +function restrict_coordinates!(v::SVector{2,Float64}, half_box_len::Float64) + return SVector( + restrict_coordinate(v[1], half_box_len), restrict_coordinate(v[2], half_box_len) + ) end -function restrict_coordinates!(p::Particle, l::Float64) - p.tmp_c = restrict_coordinates!(p.tmp_c, l) +function restrict_coordinates!(p::Particle, half_box_len::Float64) + p.tmp_c = restrict_coordinates!(p.tmp_c, half_box_len) return nothing end -function minimum_image_coordinate(value::Float64, l::Float64) - if value <= -l - value += 2 * l - elseif value > l - value -= 2 * l +function minimum_image_coordinate(value::Float64, half_box_len::Float64) + if value <= -half_box_len + value += 2 * half_box_len + elseif value > half_box_len + value -= 2 * half_box_len end return value end -function minimum_image(v::SVector{2,Float64}, l::Float64) - return SVector(minimum_image_coordinate(v[1], l), minimum_image_coordinate(v[2], l)) +function minimum_image(v::SVector{2,Float64}, half_box_len::Float64) + return SVector( + minimum_image_coordinate(v[1], half_box_len), + minimum_image_coordinate(v[2], half_box_len), + ) end function are_overlapping( - c1::SVector{2,Float64}, c2::SVector{2,Float64}, overlapping_r²::Float64, l::Float64 + c1::SVector{2,Float64}, + c2::SVector{2,Float64}, + overlapping_r²::Float64, + half_box_len::Float64, ) r⃗₁₂ = c2 - c1 # 1 -> 2 - r⃗₁₂ = minimum_image(r⃗₁₂, l) + r⃗₁₂ = minimum_image(r⃗₁₂, half_box_len) distance² = r⃗₁₂[1]^2 + r⃗₁₂[2]^2 diff --git a/src/PreVector.jl b/src/PreVector.jl index 317e3e9..72e96f5 100644 --- a/src/PreVector.jl +++ b/src/PreVector.jl @@ -4,7 +4,7 @@ mutable struct PreVector{T} last_ind::UInt64 v::Vector{T} - PreVector{T}(N) where {T} = new{T}(UInt64(0), Vector{T}(undef, N)) + PreVector{T}(n_particles) where {T} = new{T}(UInt64(0), Vector{T}(undef, n_particles)) end function push!(pv::PreVector{T}, x::T) where {T} diff --git a/src/analysis/pair_correlation_function.jl b/src/analysis/pair_correlation_function.jl index 7b0765e..b27a921 100644 --- a/src/analysis/pair_correlation_function.jl +++ b/src/analysis/pair_correlation_function.jl @@ -28,13 +28,13 @@ function pair_correlation(sol, variables) dr = (max_radius - min_radius) / n_r radius = range(min_radius, max_radius; length=n_r) - N_g = zeros(variables.N, n_r) + N_g = zeros(variables.n_particles, n_r) Threads.@threads for r_ind in 1:n_r r = radius[r_ind] - @simd for i in 1:(variables.N) - for j in 1:(variables.N) + @simd for i in 1:(variables.n_particles) + for j in 1:(variables.n_particles) if i != j for k in 1:n_last_frames frame = variables.n_snapshots - k + 1 @@ -43,7 +43,7 @@ function pair_correlation(sol, variables) r⃗₁₂ = c1 - c2 - r⃗₁₂ = minimum_image(r⃗₁₂, variables.l) + r⃗₁₂ = minimum_image(r⃗₁₂, variables.half_box_len) distance = sqrt(r⃗₁₂[1]^2 + r⃗₁₂[2]^2) @@ -62,13 +62,19 @@ function pair_correlation(sol, variables) r = radius[r_ind] tmp_g = 0.0 - @turbo for i in 1:(variables.N) + @turbo for i in 1:(variables.n_particles) tmp_g += N_g[i, r_ind] end tmp_g *= - (2 * variables.l)^2 / - ((variables.N * n_last_frames) * variables.N * 2 * π * r * dr) + (2 * variables.half_box_len)^2 / ( + (variables.n_particles * n_last_frames) * + variables.n_particles * + 2 * + π * + r * + dr + ) g[r_ind] = tmp_g end diff --git a/src/animation.jl b/src/animation.jl index e6b6aa5..6782332 100644 --- a/src/animation.jl +++ b/src/animation.jl @@ -19,7 +19,7 @@ function animate_bundle!(args) π2 = 2 * π for frame in 1:length(bundle_t) - @simd for i in 1:(args.N) + @simd for i in 1:(args.n_particles) c = bundle_c[i, frame] args.circles[][i] = Circle(Point2(c[1], c[2]), args.particle_radius) @@ -48,7 +48,7 @@ function animate_bundle!(args) println("Recording started!") else if args.debug && frame > 1 - @simd for i in 1:(args.N) + @simd for i in 1:(args.n_particles) first_ind = 2 * i - 1 second_ind = 2 * i frame_min_1 = frame - 1 @@ -96,7 +96,12 @@ function animate_after_loading(dir, animation_path, sim_consts, framerate, debug fig = Figure(; resolution=(1080, 1080)) ax = Axis( fig[1, 1]; - limits=(-sim_consts.l, sim_consts.l, -sim_consts.l, sim_consts.l), + limits=( + -sim_consts.half_box_len, + sim_consts.half_box_len, + -sim_consts.half_box_len, + sim_consts.half_box_len, + ), aspect=AxisAspect(1), xlabel=L"x", ylabel=L"y", @@ -106,11 +111,11 @@ function animate_after_loading(dir, animation_path, sim_consts, framerate, debug Colorbar(fig[1, 2]; limits=(0, 2), colormap=color_scheme, label=L"\frac{\varphi}{\pi}") - N = sim_consts.N + n_particles = sim_consts.n_particles record(fig, animation_path; framerate=framerate) do io - circles = Observable(Vector{Circle}(undef, N)) - colors = Observable(Vector{RGBAf}(undef, N)) + circles = Observable(Vector{Circle}(undef, n_particles)) + colors = Observable(Vector{RGBAf}(undef, n_particles)) segments_x = segments_y = @@ -118,14 +123,14 @@ function animate_after_loading(dir, animation_path, sim_consts, framerate, debug skin_circles = interaction_colors = skin_colors = nothing if debug - segments_x = Observable(zeros(2 * N)) - segments_y = Observable(zeros(2 * N)) + segments_x = Observable(zeros(2 * n_particles)) + segments_y = Observable(zeros(2 * n_particles)) - interaction_circles = Observable(Vector{Circle}(undef, N)) - skin_circles = Observable(Vector{Circle}(undef, N)) + interaction_circles = Observable(Vector{Circle}(undef, n_particles)) + skin_circles = Observable(Vector{Circle}(undef, n_particles)) - interaction_colors = Observable(Vector{RGBAf}(undef, N)) - skin_colors = Observable(Vector{RGBAf}(undef, N)) + interaction_colors = Observable(Vector{RGBAf}(undef, n_particles)) + skin_colors = Observable(Vector{RGBAf}(undef, n_particles)) end particle_radius = sim_consts.particle_diameter / 2 @@ -146,9 +151,9 @@ function animate_after_loading(dir, animation_path, sim_consts, framerate, debug @showprogress 1 for (n_bundle, bundle_path) in enumerate(bundle_paths) bundle::Bundle = JLD2.load_object(bundle_path) - run_vars_file = "$dir/runs/run_vars_$n_bundle.json" - if debug && isfile(run_vars_file) - skin_r::Float64 = JSON3.read(read(run_vars_file, String)).skin_r + run_params_file = "$dir/runs/run_params_$n_bundle.json" + if debug && isfile(run_params_file) + skin_r::Float64 = JSON3.read(read(run_params_file, String)).skin_r end args = (; @@ -167,7 +172,7 @@ function animate_after_loading(dir, animation_path, sim_consts, framerate, debug skin_colors, particle_radius, skin_r, - N, + n_particles, interaction_r, color_scheme, n_bundle, diff --git a/src/benchmark.jl b/src/benchmark.jl index 28390b8..8247e13 100644 --- a/src/benchmark.jl +++ b/src/benchmark.jl @@ -7,14 +7,14 @@ using ReCo function run_benchmarks( dir::String=""; - N::Int64=1000, + n_particles::Int64=1000, v::Float64=20.0, duration::Float64=2.0, n_bundle_snapshots::Int64=0, comment="", ) if length(dir) == 0 - dir = init_sim(; N=N, v=v, parent_dir="benchmark") + dir = init_sim(; n_particles=n_particles, v=v, parent_dir="benchmark") end benchmark = @benchmark run_sim( @@ -27,7 +27,7 @@ function run_benchmarks( JSON3.pretty( f, Dict( - "N" => N, + "n_particles" => n_particles, "v" => v, "duration" => duration, "n_bundle_snapshots" => n_bundle_snapshots, diff --git a/src/data.jl b/src/data.jl index 514bb21..340171f 100644 --- a/src/data.jl +++ b/src/data.jl @@ -7,11 +7,11 @@ struct Bundle φ::Matrix{Float64} end -function Bundle(N::Int64, n_snapshots::Int64) +function Bundle(n_particles::Int64, n_snapshots::Int64) return Bundle( Vector{Float64}(undef, n_snapshots), - Matrix{SVector{2,Float64}}(undef, (N, n_snapshots)), - Matrix{Float64}(undef, (N, n_snapshots)), + Matrix{SVector{2,Float64}}(undef, (n_particles, n_snapshots)), + Matrix{Float64}(undef, (n_particles, n_snapshots)), ) end @@ -33,8 +33,8 @@ function save_snapshot!( return nothing end -function set_status(dir::String, n_bundles::Int64, T::Float64) - open("$dir/status.json", "w") do f +function set_sim_state(dir::String, n_bundles::Int64, T::Float64) + open("$dir/sim_state.json", "w") do f JSON3.write(f, Dict("n_bundles" => n_bundles, "T" => round(T; digits=3))) end @@ -47,7 +47,7 @@ function save_bundle(dir::String, bundle::Bundle, n::Int64, T::Float64) JLD2.save_object("$bundles_dir/bundle_$n.jld2", bundle) - set_status(dir, n, T) + set_sim_state(dir, n, T) return nothing end \ No newline at end of file diff --git a/src/run.jl b/src/run.jl index f67b87e..a53c1ee 100644 --- a/src/run.jl +++ b/src/run.jl @@ -35,13 +35,13 @@ function run_sim( n_bundle_snapshots = min(n_snapshots, n_bundle_snapshots) - status = JSON3.read(read("$dir/status.json", String)) - n_bundles = status.n_bundles + sim_state = JSON3.read(read("$dir/sim_state.json", String)) + n_bundles = sim_state.n_bundles bundles_dir = "$dir/bundles" bundle = JLD2.load_object("$bundles_dir/bundle_$n_bundles.jld2") - particles = generate_particles(bundle, sim_consts.N) + particles = generate_particles(bundle, sim_consts.n_particles) args = ( v=sim_consts.v, @@ -56,23 +56,26 @@ function run_sim( μ=sim_consts.μ, interaction_r=sim_consts.interaction_r, interaction_r²=sim_consts.interaction_r^2, - N=sim_consts.N, - l=sim_consts.l, + n_particles=sim_consts.n_particles, + half_box_len=sim_consts.half_box_len, particle_diameter=sim_consts.particle_diameter, particles=particles, - particles_c=[particles[i].c for i in 1:(sim_consts.N)], - verlet_list=[PreVector{Int64}(sim_consts.N - 1) for i in 1:(sim_consts.N - 1)], + particles_c=[particles[i].c for i in 1:(sim_consts.n_particles)], + verlet_lists=[ + PreVector{Int64}(sim_consts.n_particles - i) for + i in 1:(sim_consts.n_particles - 1) + ], n_bundle_snapshots=n_bundle_snapshots, - bundle=Bundle(sim_consts.N, n_bundle_snapshots), - box=Box(SVector(2 * sim_consts.l, 2 * sim_consts.l), skin_r), + bundle=Bundle(sim_consts.n_particles, n_bundle_snapshots), + box=Box(SVector(2 * sim_consts.half_box_len, 2 * sim_consts.half_box_len), skin_r), ) - T0::Float64 = status.T + T0::Float64 = sim_state.T T = T0 + duration start_datetime = now() - run_vars = (; + run_params = (; # Input duration, snapshot_at, @@ -98,8 +101,8 @@ function run_sim( if n_bundle_snapshots > 0 save_data = true - open("$runs_dir/run_vars_$next_bundle.json", "w") do f - JSON3.write(f, run_vars) + open("$runs_dir/run_params_$next_bundle.json", "w") do f + JSON3.write(f, run_params) end else save_data = false diff --git a/src/setup.jl b/src/setup.jl index a0f4d65..376b064 100644 --- a/src/setup.jl +++ b/src/setup.jl @@ -2,12 +2,14 @@ using Distributions: Uniform using Dates: now using JSON3: JSON3 -function initial_particle_grid_pos(i::Int64, j::Int64, grid_box_width::Float64, l::Float64) - term = -0.5 * grid_box_width - l +function initial_particle_grid_pos( + i::Int64, j::Int64, grid_box_width::Float64, half_box_len::Float64 +) + term = -0.5 * grid_box_width - half_box_len return SVector(i * grid_box_width + term, j * grid_box_width + term) end -function generate_particles(grid_n::Int64, grid_box_width::Float64, l::Float64) +function generate_particles(grid_n::Int64, grid_box_width::Float64, half_box_len::Float64) particles = Vector{Particle}(undef, grid_n^2) id = 1 @@ -15,7 +17,9 @@ function generate_particles(grid_n::Int64, grid_box_width::Float64, l::Float64) for i in 1:grid_n for j in 1:grid_n particles[id] = Particle( - id, initial_particle_grid_pos(i, j, grid_box_width, l), rand(Uniform(-π, π)) + id, + initial_particle_grid_pos(i, j, grid_box_width, half_box_len), + rand(Uniform(-π, π)), ) id += 1 @@ -25,10 +29,10 @@ function generate_particles(grid_n::Int64, grid_box_width::Float64, l::Float64) return particles end -function generate_particles(bundle::Bundle, N::Int64) - particles = Vector{Particle}(undef, N) +function generate_particles(bundle::Bundle, n_particles::Int64) + particles = Vector{Particle}(undef, n_particles) - @simd for id in 1:N + @simd for id in 1:n_particles particles[id] = Particle(id, bundle.c[id, end], bundle.φ[id, end]) end @@ -36,14 +40,14 @@ function generate_particles(bundle::Bundle, N::Int64) end function init_sim(; - N::Int64, + n_particles::Int64, v::Float64, δt::Float64=1e-5, packing_ratio::Float64=0.5, comment::String="", parent_dir::String="", ) - @assert N > 0 + @assert n_particles > 0 @assert v >= 0 @assert δt in 1e-7:1e-7:1e-4 @assert packing_ratio > 0 @@ -56,16 +60,16 @@ function init_sim(; ϵ = 100.0 interaction_r = 2^(1 / 6) * σ - grid_n = round(Int64, ceil(sqrt(N))) + grid_n = round(Int64, ceil(sqrt(n_particles))) - N = grid_n^2 + n_particles = grid_n^2 - l = sqrt(N * π / packing_ratio) * σ / 4 - grid_box_width = 2 * l / grid_n + half_box_len = sqrt(n_particles * π / packing_ratio) * σ / 4 + grid_box_width = 2 * half_box_len / grid_n sim_consts = (; # Input - N, + n_particles, v, δt, packing_ratio, @@ -78,12 +82,12 @@ function init_sim(; ϵ, interaction_r, grid_n, - l, + half_box_len, grid_box_width, ) - particles = generate_particles(grid_n, grid_box_width, l) - bundle = Bundle(N, 1) + particles = generate_particles(grid_n, grid_box_width, half_box_len) + bundle = Bundle(n_particles, 1) save_snapshot!(bundle, 1, 0.0, particles) particles = nothing @@ -92,7 +96,7 @@ function init_sim(; if length(parent_dir) > 0 dir *= "/$parent_dir" end - dir *= "/$(now())_N=$(N)_v=$(v)_#$(rand(1000:9999))" + dir *= "/$(now())_N=$(n_particles)_v=$(v)_#$(rand(1000:9999))" if length(comment) > 0 dir *= "_$comment" end diff --git a/src/shape.jl b/src/shape.jl index 5019125..66e5f39 100644 --- a/src/shape.jl +++ b/src/shape.jl @@ -1,30 +1,30 @@ using StaticArrays: SVector, SMatrix using LinearAlgebra: eigvals, Hermitian -function project_to_unit_circle(x::Float64, l::Float64) - φ = (x + l) * π / l +function project_to_unit_circle(x::Float64, half_box_len::Float64) + φ = (x + half_box_len) * π / half_box_len si, co = sincos(φ) return SVector(co, si) end -function project_back_from_unit_circle(θ::T, l::Float64) where {T<:Real} - x = θ * l / π - l - return restrict_coordinate(x, l) +function project_back_from_unit_circle(θ::T, half_box_len::Float64) where {T<:Real} + x = θ * half_box_len / π - half_box_len + return restrict_coordinate(x, half_box_len) end -function center_of_mass(particles::Vector{Particle}, l::Float64) +function center_of_mass(particles::Vector{Particle}, half_box_len::Float64) x_proj_sum = SVector(0.0, 0.0) y_proj_sum = SVector(0.0, 0.0) for p in particles - x_proj_sum += project_to_unit_circle(p.c[1], l) - y_proj_sum += project_to_unit_circle(p.c[2], l) + x_proj_sum += project_to_unit_circle(p.c[1], half_box_len) + y_proj_sum += project_to_unit_circle(p.c[2], half_box_len) end # Prevent for example atan(1e-16, 1e-15) != 0 with rounding digits = 5 - # No need for 1/N with atan + # No need for 1/n_particles with atan # If proj is (0, 0) then COM is 0 or L or -L. Here, 0 is choosen with θ = π if round(x_proj_sum[1]; digits=digits) == round(x_proj_sum[2]; digits=digits) == 0 x_θ = π @@ -38,21 +38,21 @@ function center_of_mass(particles::Vector{Particle}, l::Float64) y_θ = atan(y_proj_sum[2], y_proj_sum[1]) end - COM_x = project_back_from_unit_circle(x_θ, l) - COM_y = project_back_from_unit_circle(y_θ, l) + COM_x = project_back_from_unit_circle(x_θ, half_box_len) + COM_y = project_back_from_unit_circle(y_θ, half_box_len) return SVector(COM_x, COM_y) end -function gyration_tensor(particles::Vector{Particle}, l::Float64) - COM = center_of_mass(particles, l) +function gyration_tensor(particles::Vector{Particle}, half_box_len::Float64) + COM = center_of_mass(particles, half_box_len) S11 = 0.0 S12 = 0.0 S22 = 0.0 for p in particles - shifted_c = restrict_coordinates!(p.c - COM, l) + shifted_c = restrict_coordinates!(p.c - COM, half_box_len) S11 += shifted_c[1]^2 S12 += shifted_c[1] * shifted_c[2] @@ -62,7 +62,7 @@ function gyration_tensor(particles::Vector{Particle}, l::Float64) return Hermitian(SMatrix{2,2}(S11, S12, S12, S22)) end -function gyration_tensor_eigvals_ratio(particles::Vector{Particle}, l::Float64) - ev = eigvals(gyration_tensor(particles, l)) # Eigenvalues are sorted +function gyration_tensor_eigvals_ratio(particles::Vector{Particle}, half_box_len::Float64) + ev = eigvals(gyration_tensor(particles, half_box_len)) # Eigenvalues are sorted return ev[1] / ev[2] end \ No newline at end of file diff --git a/src/simulation.jl b/src/simulation.jl index add6857..eb362e4 100644 --- a/src/simulation.jl +++ b/src/simulation.jl @@ -8,29 +8,29 @@ import Base: wait rand_normal01() = rand(Normal(0, 1)) -function push_to_verlet_list(verlet_list, i, j) +function push_to_verlet_list!(verlet_lists, i, j) if i < j - push!(verlet_list[i], Int64(j)) + push!(verlet_lists[i], Int64(j)) else - push!(verlet_list[j], Int64(i)) + push!(verlet_lists[j], Int64(i)) end return nothing end -function update_verlet_list!(args, cl) - @simd for pv in args.verlet_list - reset!(pv) +function update_verlet_lists!(args, cl) + @simd for pre_vec in args.verlet_lists + reset!(pre_vec) end - @simd for i in 1:(args.N) + @simd for i in 1:(args.n_particles) args.particles_c[i] = args.particles[i].c end cl = UpdateCellList!(args.particles_c, args.box, cl; parallel=false) map_pairwise!( - (x, y, i, j, d2, output) -> push_to_verlet_list(args.verlet_list, i, j), + (x, y, i, j, d2, output) -> push_to_verlet_list!(args.verlet_lists, i, j), nothing, args.box, cl; @@ -41,16 +41,16 @@ function update_verlet_list!(args, cl) end function euler!(args) - for id1 in 1:(args.N - 1) + for id1 in 1:(args.n_particles - 1) p1 = args.particles[id1] p1_c = p1.c - verlet_list = args.verlet_list[id1] + verlet_list = args.verlet_lists[id1] for id2 in view(verlet_list.v, 1:(verlet_list.last_ind)) p2 = args.particles[id2] overlapping, r⃗₁₂, distance² = are_overlapping( - p1_c, p2.c, args.interaction_r², args.l + p1_c, p2.c, args.interaction_r², args.half_box_len ) if overlapping @@ -72,7 +72,7 @@ function euler!(args) p.φ += args.c₄ * rand_normal01() - restrict_coordinates!(p, args.l) + restrict_coordinates!(p, args.half_box_len) p.c = p.tmp_c end @@ -98,7 +98,7 @@ function simulate( task::Union{Task,Nothing} = nothing cl = CellList(args.particles_c, args.box; parallel=false) - cl = update_verlet_list!(args, cl) + cl = update_verlet_lists!(args, cl) start_time = now() println("Started simulation at $start_time.") @@ -121,7 +121,7 @@ function simulate( end if integration_step % n_steps_before_verlet_list_update == 0 - cl = update_verlet_list!(args, cl) + cl = update_verlet_lists!(args, cl) end euler!(args)