mo8it/ReCo.jl
1
0
Fork 0
mirror of https://gitlab.rlp.net/mobitar/ReCo.jl.git synced 2024-12-21 00:51:21 +00:00
Code Issues Releases Activity

Renaming and verlet list shrink

Browse source
This commit is contained in:
MoBit 2021-12-10 03:16:18 +01:00
parent 5001c315ce
commit 50feaee469
11 changed files with 138 additions and 113 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

1
Project.toml
View file

@ -22,7 +22,6 @@ ProfileView = "c46f51b8-102a-5cf2-8d2c-8597cb0e0da7"
ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca" ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca"
Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
ReinforcementLearning = "158674fc-8238-5cab-b5ba-03dfc80d1318" ReinforcementLearning = "158674fc-8238-5cab-b5ba-03dfc80d1318"
Revise = "295af30f-e4ad-537b-8983-00126c2a3abe"
StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

44
src/Particle.jl
View file

@ -13,46 +13,54 @@ mutable struct Particle
end end
end end
function restrict_coordinate(value::Float64, l::Float64) function restrict_coordinate(value::Float64, half_box_len::Float64)
if value < -l if value < -half_box_len
value += 2 * l value += 2 * half_box_len
elseif value >= l elseif value >= half_box_len
value -= 2 * l value -= 2 * half_box_len
end end
return value return value
end end
function restrict_coordinates!(v::SVector{2,Float64}, l::Float64) function restrict_coordinates!(v::SVector{2,Float64}, half_box_len::Float64)
return SVector(restrict_coordinate(v[1], l), restrict_coordinate(v[2], l)) return SVector(
restrict_coordinate(v[1], half_box_len), restrict_coordinate(v[2], half_box_len)
)
end end
function restrict_coordinates!(p::Particle, l::Float64) function restrict_coordinates!(p::Particle, half_box_len::Float64)
p.tmp_c = restrict_coordinates!(p.tmp_c, l) p.tmp_c = restrict_coordinates!(p.tmp_c, half_box_len)
return nothing return nothing
end end
function minimum_image_coordinate(value::Float64, l::Float64) function minimum_image_coordinate(value::Float64, half_box_len::Float64)
if value <= -l if value <= -half_box_len
value += 2 * l value += 2 * half_box_len
elseif value > l elseif value > half_box_len
value -= 2 * l value -= 2 * half_box_len
end end
return value return value
end end
function minimum_image(v::SVector{2,Float64}, l::Float64) function minimum_image(v::SVector{2,Float64}, half_box_len::Float64)
return SVector(minimum_image_coordinate(v[1], l), minimum_image_coordinate(v[2], l)) return SVector(
minimum_image_coordinate(v[1], half_box_len),
minimum_image_coordinate(v[2], half_box_len),
)
end end
function are_overlapping( 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⃗₁₂ = c2 - c1 # 1 -> 2
r⃗₁₂ = minimum_image(r⃗₁₂, l) r⃗₁₂ = minimum_image(r⃗₁₂, half_box_len)
distance² = r⃗₁₂[1]^2 + r⃗₁₂[2]^2 distance² = r⃗₁₂[1]^2 + r⃗₁₂[2]^2

2
src/PreVector.jl
View file

@ -4,7 +4,7 @@ mutable struct PreVector{T}
last_ind::UInt64 last_ind::UInt64
v::Vector{T} 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 end
function push!(pv::PreVector{T}, x::T) where {T} function push!(pv::PreVector{T}, x::T) where {T}

20
src/analysis/pair_correlation_function.jl
View file

@ -28,13 +28,13 @@ function pair_correlation(sol, variables)
dr = (max_radius - min_radius) / n_r dr = (max_radius - min_radius) / n_r
radius = range(min_radius, max_radius; length=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 Threads.@threads for r_ind in 1:n_r
r = radius[r_ind] r = radius[r_ind]
@simd for i in 1:(variables.N) @simd for i in 1:(variables.n_particles)
for j in 1:(variables.N) for j in 1:(variables.n_particles)
if i != j if i != j
for k in 1:n_last_frames for k in 1:n_last_frames
frame = variables.n_snapshots - k + 1 frame = variables.n_snapshots - k + 1
@ -43,7 +43,7 @@ function pair_correlation(sol, variables)
r⃗₁₂ = c1 - c2 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) distance = sqrt(r⃗₁₂[1]^2 + r⃗₁₂[2]^2)
@ -62,13 +62,19 @@ function pair_correlation(sol, variables)
r = radius[r_ind] r = radius[r_ind]
tmp_g = 0.0 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] tmp_g += N_g[i, r_ind]
end end
tmp_g *= tmp_g *=
(2 * variables.l)^2 / (2 * variables.half_box_len)^2 / (
((variables.N * n_last_frames) * variables.N * 2 * π * r * dr) (variables.n_particles * n_last_frames) *
variables.n_particles *
2 *
π *
r *
dr
)
g[r_ind] = tmp_g g[r_ind] = tmp_g
end end

37
src/animation.jl
View file

@ -19,7 +19,7 @@ function animate_bundle!(args)
π2 = 2 * π π2 = 2 * π
for frame in 1:length(bundle_t) 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] c = bundle_c[i, frame]
args.circles[][i] = Circle(Point2(c[1], c[2]), args.particle_radius) args.circles[][i] = Circle(Point2(c[1], c[2]), args.particle_radius)
@ -48,7 +48,7 @@ function animate_bundle!(args)
println("Recording started!") println("Recording started!")
else else
if args.debug && frame > 1 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 first_ind = 2 * i - 1
second_ind = 2 * i second_ind = 2 * i
frame_min_1 = frame - 1 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)) fig = Figure(; resolution=(1080, 1080))
ax = Axis( ax = Axis(
fig[1, 1]; 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), aspect=AxisAspect(1),
xlabel=L"x", xlabel=L"x",
ylabel=L"y", 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}") 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 record(fig, animation_path; framerate=framerate) do io
circles = Observable(Vector{Circle}(undef, N)) circles = Observable(Vector{Circle}(undef, n_particles))
colors = Observable(Vector{RGBAf}(undef, N)) colors = Observable(Vector{RGBAf}(undef, n_particles))
segments_x = segments_x =
segments_y = segments_y =
@ -118,14 +123,14 @@ function animate_after_loading(dir, animation_path, sim_consts, framerate, debug
skin_circles = interaction_colors = skin_colors = nothing skin_circles = interaction_colors = skin_colors = nothing
if debug if debug
segments_x = Observable(zeros(2 * N)) segments_x = Observable(zeros(2 * n_particles))
segments_y = Observable(zeros(2 * N)) segments_y = Observable(zeros(2 * n_particles))
interaction_circles = Observable(Vector{Circle}(undef, N)) interaction_circles = Observable(Vector{Circle}(undef, n_particles))
skin_circles = Observable(Vector{Circle}(undef, N)) skin_circles = Observable(Vector{Circle}(undef, n_particles))
interaction_colors = Observable(Vector{RGBAf}(undef, N)) interaction_colors = Observable(Vector{RGBAf}(undef, n_particles))
skin_colors = Observable(Vector{RGBAf}(undef, N)) skin_colors = Observable(Vector{RGBAf}(undef, n_particles))
end end
particle_radius = sim_consts.particle_diameter / 2 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) @showprogress 1 for (n_bundle, bundle_path) in enumerate(bundle_paths)
bundle::Bundle = JLD2.load_object(bundle_path) bundle::Bundle = JLD2.load_object(bundle_path)
run_vars_file = "$dir/runs/run_vars_$n_bundle.json" run_params_file = "$dir/runs/run_params_$n_bundle.json"
if debug && isfile(run_vars_file) if debug && isfile(run_params_file)
skin_r::Float64 = JSON3.read(read(run_vars_file, String)).skin_r skin_r::Float64 = JSON3.read(read(run_params_file, String)).skin_r
end end
args = (; args = (;
@ -167,7 +172,7 @@ function animate_after_loading(dir, animation_path, sim_consts, framerate, debug
skin_colors, skin_colors,
particle_radius, particle_radius,
skin_r, skin_r,
N, n_particles,
interaction_r, interaction_r,
color_scheme, color_scheme,
n_bundle, n_bundle,

6
src/benchmark.jl
View file

@ -7,14 +7,14 @@ using ReCo
function run_benchmarks( function run_benchmarks(
dir::String=""; dir::String="";
N::Int64=1000, n_particles::Int64=1000,
v::Float64=20.0, v::Float64=20.0,
duration::Float64=2.0, duration::Float64=2.0,
n_bundle_snapshots::Int64=0, n_bundle_snapshots::Int64=0,
comment="", comment="",
) )
if length(dir) == 0 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 end
benchmark = @benchmark run_sim( benchmark = @benchmark run_sim(
@ -27,7 +27,7 @@ function run_benchmarks(
JSON3.pretty( JSON3.pretty(
f, f,
Dict( Dict(
"N" => N, "n_particles" => n_particles,
"v" => v, "v" => v,
"duration" => duration, "duration" => duration,
"n_bundle_snapshots" => n_bundle_snapshots, "n_bundle_snapshots" => n_bundle_snapshots,

12
src/data.jl
View file

@ -7,11 +7,11 @@ struct Bundle
φ::Matrix{Float64} φ::Matrix{Float64}
end end
function Bundle(N::Int64, n_snapshots::Int64) function Bundle(n_particles::Int64, n_snapshots::Int64)
return Bundle( return Bundle(
Vector{Float64}(undef, n_snapshots), Vector{Float64}(undef, n_snapshots),
Matrix{SVector{2,Float64}}(undef, (N, n_snapshots)), Matrix{SVector{2,Float64}}(undef, (n_particles, n_snapshots)),
Matrix{Float64}(undef, (N, n_snapshots)), Matrix{Float64}(undef, (n_particles, n_snapshots)),
) )
end end
@ -33,8 +33,8 @@ function save_snapshot!(
return nothing return nothing
end end
function set_status(dir::String, n_bundles::Int64, T::Float64) function set_sim_state(dir::String, n_bundles::Int64, T::Float64)
open("$dir/status.json", "w") do f open("$dir/sim_state.json", "w") do f
JSON3.write(f, Dict("n_bundles" => n_bundles, "T" => round(T; digits=3))) JSON3.write(f, Dict("n_bundles" => n_bundles, "T" => round(T; digits=3)))
end 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) JLD2.save_object("$bundles_dir/bundle_$n.jld2", bundle)
set_status(dir, n, T) set_sim_state(dir, n, T)
return nothing return nothing
end end

29
src/run.jl
View file

@ -35,13 +35,13 @@ function run_sim(
n_bundle_snapshots = min(n_snapshots, n_bundle_snapshots) n_bundle_snapshots = min(n_snapshots, n_bundle_snapshots)
status = JSON3.read(read("$dir/status.json", String)) sim_state = JSON3.read(read("$dir/sim_state.json", String))
n_bundles = status.n_bundles n_bundles = sim_state.n_bundles
bundles_dir = "$dir/bundles" bundles_dir = "$dir/bundles"
bundle = JLD2.load_object("$bundles_dir/bundle_$n_bundles.jld2") 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 = ( args = (
v=sim_consts.v, v=sim_consts.v,
@ -56,23 +56,26 @@ function run_sim(
μ=sim_consts.μ, μ=sim_consts.μ,
interaction_r=sim_consts.interaction_r, interaction_r=sim_consts.interaction_r,
interaction_r²=sim_consts.interaction_r^2, interaction_r²=sim_consts.interaction_r^2,
N=sim_consts.N, n_particles=sim_consts.n_particles,
l=sim_consts.l, half_box_len=sim_consts.half_box_len,
particle_diameter=sim_consts.particle_diameter, particle_diameter=sim_consts.particle_diameter,
particles=particles, particles=particles,
particles_c=[particles[i].c for i in 1:(sim_consts.N)], particles_c=[particles[i].c for i in 1:(sim_consts.n_particles)],
verlet_list=[PreVector{Int64}(sim_consts.N - 1) for i in 1:(sim_consts.N - 1)], verlet_lists=[
PreVector{Int64}(sim_consts.n_particles - i) for
i in 1:(sim_consts.n_particles - 1)
],
n_bundle_snapshots=n_bundle_snapshots, n_bundle_snapshots=n_bundle_snapshots,
bundle=Bundle(sim_consts.N, n_bundle_snapshots), bundle=Bundle(sim_consts.n_particles, n_bundle_snapshots),
box=Box(SVector(2 * sim_consts.l, 2 * sim_consts.l), skin_r), 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 T = T0 + duration
start_datetime = now() start_datetime = now()
run_vars = (; run_params = (;
# Input # Input
duration, duration,
snapshot_at, snapshot_at,
@ -98,8 +101,8 @@ function run_sim(
if n_bundle_snapshots > 0 if n_bundle_snapshots > 0
save_data = true save_data = true
open("$runs_dir/run_vars_$next_bundle.json", "w") do f open("$runs_dir/run_params_$next_bundle.json", "w") do f
JSON3.write(f, run_vars) JSON3.write(f, run_params)
end end
else else
save_data = false save_data = false

40
src/setup.jl
View file

@ -2,12 +2,14 @@ using Distributions: Uniform
using Dates: now using Dates: now
using JSON3: JSON3 using JSON3: JSON3
function initial_particle_grid_pos(i::Int64, j::Int64, grid_box_width::Float64, l::Float64) function initial_particle_grid_pos(
term = -0.5 * grid_box_width - l 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) return SVector(i * grid_box_width + term, j * grid_box_width + term)
end 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) particles = Vector{Particle}(undef, grid_n^2)
id = 1 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 i in 1:grid_n
for j in 1:grid_n for j in 1:grid_n
particles[id] = Particle( 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 id += 1
@ -25,10 +29,10 @@ function generate_particles(grid_n::Int64, grid_box_width::Float64, l::Float64)
return particles return particles
end end
function generate_particles(bundle::Bundle, N::Int64) function generate_particles(bundle::Bundle, n_particles::Int64)
particles = Vector{Particle}(undef, N) 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]) particles[id] = Particle(id, bundle.c[id, end], bundle.φ[id, end])
end end
@ -36,14 +40,14 @@ function generate_particles(bundle::Bundle, N::Int64)
end end
function init_sim(; function init_sim(;
N::Int64, n_particles::Int64,
v::Float64, v::Float64,
δt::Float64=1e-5, δt::Float64=1e-5,
packing_ratio::Float64=0.5, packing_ratio::Float64=0.5,
comment::String="", comment::String="",
parent_dir::String="", parent_dir::String="",
) )
@assert N > 0 @assert n_particles > 0
@assert v >= 0 @assert v >= 0
@assert δt in 1e-7:1e-7:1e-4 @assert δt in 1e-7:1e-7:1e-4
@assert packing_ratio > 0 @assert packing_ratio > 0
@ -56,16 +60,16 @@ function init_sim(;
ϵ = 100.0 ϵ = 100.0
interaction_r = 2^(1 / 6) * σ 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 half_box_len = sqrt(n_particles * π / packing_ratio) * σ / 4
grid_box_width = 2 * l / grid_n grid_box_width = 2 * half_box_len / grid_n
sim_consts = (; sim_consts = (;
# Input # Input
N, n_particles,
v, v,
δt, δt,
packing_ratio, packing_ratio,
@ -78,12 +82,12 @@ function init_sim(;
ϵ, ϵ,
interaction_r, interaction_r,
grid_n, grid_n,
l, half_box_len,
grid_box_width, grid_box_width,
) )
particles = generate_particles(grid_n, grid_box_width, l) particles = generate_particles(grid_n, grid_box_width, half_box_len)
bundle = Bundle(N, 1) bundle = Bundle(n_particles, 1)
save_snapshot!(bundle, 1, 0.0, particles) save_snapshot!(bundle, 1, 0.0, particles)
particles = nothing particles = nothing
@ -92,7 +96,7 @@ function init_sim(;
if length(parent_dir) > 0 if length(parent_dir) > 0
dir *= "/$parent_dir" dir *= "/$parent_dir"
end end
dir *= "/$(now())_N=$(N)_v=$(v)_#$(rand(1000:9999))" dir *= "/$(now())_N=$(n_particles)_v=$(v)_#$(rand(1000:9999))"
if length(comment) > 0 if length(comment) > 0
dir *= "_$comment" dir *= "_$comment"
end end

32
src/shape.jl
View file

@ -1,30 +1,30 @@
using StaticArrays: SVector, SMatrix using StaticArrays: SVector, SMatrix
using LinearAlgebra: eigvals, Hermitian using LinearAlgebra: eigvals, Hermitian
function project_to_unit_circle(x::Float64, l::Float64) function project_to_unit_circle(x::Float64, half_box_len::Float64)
φ = (x + l) * π / l φ = (x + half_box_len) * π / half_box_len
si, co = sincos(φ) si, co = sincos(φ)
return SVector(co, si) return SVector(co, si)
end end
function project_back_from_unit_circle(θ::T, l::Float64) where {T<:Real} function project_back_from_unit_circle(θ::T, half_box_len::Float64) where {T<:Real}
x = θ * l / π - l x = θ * half_box_len / π - half_box_len
return restrict_coordinate(x, l) return restrict_coordinate(x, half_box_len)
end 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) x_proj_sum = SVector(0.0, 0.0)
y_proj_sum = SVector(0.0, 0.0) y_proj_sum = SVector(0.0, 0.0)
for p in particles for p in particles
x_proj_sum += project_to_unit_circle(p.c[1], l) x_proj_sum += project_to_unit_circle(p.c[1], half_box_len)
y_proj_sum += project_to_unit_circle(p.c[2], l) y_proj_sum += project_to_unit_circle(p.c[2], half_box_len)
end end
# Prevent for example atan(1e-16, 1e-15) != 0 with rounding # Prevent for example atan(1e-16, 1e-15) != 0 with rounding
digits = 5 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 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 if round(x_proj_sum[1]; digits=digits) == round(x_proj_sum[2]; digits=digits) == 0
x_θ = π x_θ = π
@ -38,21 +38,21 @@ function center_of_mass(particles::Vector{Particle}, l::Float64)
y_θ = atan(y_proj_sum[2], y_proj_sum[1]) y_θ = atan(y_proj_sum[2], y_proj_sum[1])
end end
COM_x = project_back_from_unit_circle(x_θ, l) COM_x = project_back_from_unit_circle(x_θ, half_box_len)
COM_y = project_back_from_unit_circle(y_θ, l) COM_y = project_back_from_unit_circle(y_θ, half_box_len)
return SVector(COM_x, COM_y) return SVector(COM_x, COM_y)
end end
function gyration_tensor(particles::Vector{Particle}, l::Float64) function gyration_tensor(particles::Vector{Particle}, half_box_len::Float64)
COM = center_of_mass(particles, l) COM = center_of_mass(particles, half_box_len)
S11 = 0.0 S11 = 0.0
S12 = 0.0 S12 = 0.0
S22 = 0.0 S22 = 0.0
for p in particles 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 S11 += shifted_c[1]^2
S12 += shifted_c[1] * shifted_c[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)) return Hermitian(SMatrix{2,2}(S11, S12, S12, S22))
end end
function gyration_tensor_eigvals_ratio(particles::Vector{Particle}, l::Float64) function gyration_tensor_eigvals_ratio(particles::Vector{Particle}, half_box_len::Float64)
ev = eigvals(gyration_tensor(particles, l)) # Eigenvalues are sorted ev = eigvals(gyration_tensor(particles, half_box_len)) # Eigenvalues are sorted
return ev[1] / ev[2] return ev[1] / ev[2]
end end

28
src/simulation.jl
View file

@ -8,29 +8,29 @@ import Base: wait
rand_normal01() = rand(Normal(0, 1)) 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 if i < j
push!(verlet_list[i], Int64(j)) push!(verlet_lists[i], Int64(j))
else else
push!(verlet_list[j], Int64(i)) push!(verlet_lists[j], Int64(i))
end end
return nothing return nothing
end end
function update_verlet_list!(args, cl) function update_verlet_lists!(args, cl)
@simd for pv in args.verlet_list @simd for pre_vec in args.verlet_lists
reset!(pv) reset!(pre_vec)
end end
@simd for i in 1:(args.N) @simd for i in 1:(args.n_particles)
args.particles_c[i] = args.particles[i].c args.particles_c[i] = args.particles[i].c
end end
cl = UpdateCellList!(args.particles_c, args.box, cl; parallel=false) cl = UpdateCellList!(args.particles_c, args.box, cl; parallel=false)
map_pairwise!( 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, nothing,
args.box, args.box,
cl; cl;
@ -41,16 +41,16 @@ function update_verlet_list!(args, cl)
end end
function euler!(args) function euler!(args)
for id1 in 1:(args.N - 1) for id1 in 1:(args.n_particles - 1)
p1 = args.particles[id1] p1 = args.particles[id1]
p1_c = p1.c 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)) for id2 in view(verlet_list.v, 1:(verlet_list.last_ind))
p2 = args.particles[id2] p2 = args.particles[id2]
overlapping, r⃗₁₂, distance² = are_overlapping( 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 if overlapping
@ -72,7 +72,7 @@ function euler!(args)
p.φ += args.c₄ * rand_normal01() p.φ += args.c₄ * rand_normal01()
restrict_coordinates!(p, args.l) restrict_coordinates!(p, args.half_box_len)
p.c = p.tmp_c p.c = p.tmp_c
end end
@ -98,7 +98,7 @@ function simulate(
task::Union{Task,Nothing} = nothing task::Union{Task,Nothing} = nothing
cl = CellList(args.particles_c, args.box; parallel=false) cl = CellList(args.particles_c, args.box; parallel=false)
cl = update_verlet_list!(args, cl) cl = update_verlet_lists!(args, cl)
start_time = now() start_time = now()
println("Started simulation at $start_time.") println("Started simulation at $start_time.")
@ -121,7 +121,7 @@ function simulate(
end end
if integration_step % n_steps_before_verlet_list_update == 0 if integration_step % n_steps_before_verlet_list_update == 0
cl = update_verlet_list!(args, cl) cl = update_verlet_lists!(args, cl)
end end
euler!(args) euler!(args)