2021-12-07 03:19:35 +00:00
|
|
|
using Dates: Dates, now
|
2021-12-02 22:13:54 +00:00
|
|
|
using ProgressMeter: @showprogress
|
|
|
|
using LoopVectorization: @turbo
|
|
|
|
using Distributions: Normal
|
2021-12-07 03:19:35 +00:00
|
|
|
using CellListMap: Box, CellList, map_pairwise!, UpdateCellList!
|
|
|
|
|
2021-12-02 22:13:54 +00:00
|
|
|
import Base: wait
|
|
|
|
|
2021-11-10 14:41:04 +00:00
|
|
|
rand_normal01() = rand(Normal(0, 1))
|
|
|
|
|
2021-12-07 03:19:35 +00:00
|
|
|
function push_to_verlet_list(verlet_list, i, j)
|
|
|
|
if i < j
|
|
|
|
push!(verlet_list[i], Int64(j))
|
|
|
|
else
|
|
|
|
push!(verlet_list[j], Int64(i))
|
|
|
|
end
|
|
|
|
|
|
|
|
return nothing
|
|
|
|
end
|
|
|
|
|
|
|
|
function update_verlet_list!(args, cl)
|
2021-11-10 14:41:04 +00:00
|
|
|
@simd for pv in args.verlet_list
|
|
|
|
reset!(pv)
|
|
|
|
end
|
2021-11-16 21:26:08 +00:00
|
|
|
|
2021-12-07 03:19:35 +00:00
|
|
|
for i in 1:(args.N)
|
|
|
|
args.particles_c[i] = SVector(args.particles[i].c[1], args.particles[i].c[2])
|
|
|
|
end
|
2021-11-10 14:41:04 +00:00
|
|
|
|
2021-12-07 03:19:35 +00:00
|
|
|
cl = UpdateCellList!(args.particles_c, args.box, cl; parallel=false)
|
2021-11-10 14:41:04 +00:00
|
|
|
|
2021-12-07 03:19:35 +00:00
|
|
|
map_pairwise!(
|
|
|
|
(x, y, i, j, d2, output) -> push_to_verlet_list(args.verlet_list, i, j),
|
|
|
|
nothing,
|
|
|
|
args.box,
|
|
|
|
cl;
|
|
|
|
parallel=false,
|
|
|
|
)
|
2021-11-16 21:26:08 +00:00
|
|
|
|
2021-12-07 03:19:35 +00:00
|
|
|
return cl
|
2021-11-10 14:41:04 +00:00
|
|
|
end
|
|
|
|
|
|
|
|
function euler!(args)
|
2021-11-11 00:25:06 +00:00
|
|
|
for i in 1:(args.N - 1)
|
2021-12-07 00:52:34 +00:00
|
|
|
p1 = args.particles[i]
|
|
|
|
verlet_list = args.verlet_list[p1.id]
|
2021-11-10 14:41:04 +00:00
|
|
|
|
2021-11-16 21:26:08 +00:00
|
|
|
for j in 1:(verlet_list.last_ind)
|
2021-11-11 00:25:06 +00:00
|
|
|
p2 = args.particles[verlet_list.v[j]]
|
2021-11-10 14:41:04 +00:00
|
|
|
|
2021-11-16 21:26:08 +00:00
|
|
|
overlapping, r⃗₁₂, distance² = are_overlapping(
|
2021-12-07 00:52:34 +00:00
|
|
|
p1, p2, args.interaction_r², args.l
|
2021-11-16 21:26:08 +00:00
|
|
|
)
|
2021-11-10 14:41:04 +00:00
|
|
|
|
|
|
|
if overlapping
|
|
|
|
c = args.c₁ / (distance²^4) * (args.c₂ / (distance²^3) - 1)
|
2021-11-16 21:26:08 +00:00
|
|
|
|
2021-11-19 00:21:41 +00:00
|
|
|
@turbo for k in 1:2
|
2021-11-26 02:35:39 +00:00
|
|
|
dc = c * r⃗₁₂[k]
|
2021-11-16 21:26:08 +00:00
|
|
|
|
2021-12-07 00:52:34 +00:00
|
|
|
p1.tmp_c[k] -= dc
|
2021-11-26 02:35:39 +00:00
|
|
|
p2.tmp_c[k] += dc
|
2021-11-10 14:41:04 +00:00
|
|
|
end
|
|
|
|
end
|
|
|
|
end
|
2021-11-11 00:25:06 +00:00
|
|
|
end
|
2021-11-10 14:41:04 +00:00
|
|
|
|
2021-11-11 00:25:06 +00:00
|
|
|
@simd for p in args.particles
|
2021-11-10 14:41:04 +00:00
|
|
|
e = SVector(cos(p.φ), sin(p.φ))
|
|
|
|
|
2021-11-19 00:21:41 +00:00
|
|
|
@turbo for i in 1:2
|
2021-11-10 14:41:04 +00:00
|
|
|
p.tmp_c[i] += args.vδt * e[i] + args.c₃ * rand_normal01()
|
|
|
|
end
|
|
|
|
|
|
|
|
p.φ += args.c₄ * rand_normal01()
|
|
|
|
|
2021-12-07 00:36:47 +00:00
|
|
|
restrict_coordinates!(p, args.l)
|
2021-11-10 14:41:04 +00:00
|
|
|
|
2021-11-11 00:25:06 +00:00
|
|
|
update!(p)
|
2021-11-10 14:41:04 +00:00
|
|
|
end
|
|
|
|
|
|
|
|
return nothing
|
|
|
|
end
|
|
|
|
|
2021-11-26 13:50:47 +00:00
|
|
|
wait(n::Nothing) = n
|
|
|
|
|
2021-11-16 21:26:08 +00:00
|
|
|
function simulate(
|
|
|
|
args,
|
2021-11-10 23:24:07 +00:00
|
|
|
δt::Float64,
|
2021-11-26 05:55:27 +00:00
|
|
|
T0::Float64,
|
|
|
|
T::Float64,
|
2021-11-10 23:24:07 +00:00
|
|
|
n_steps_before_verlet_list_update::Int64,
|
2021-11-26 05:55:27 +00:00
|
|
|
n_steps_before_snapshot::Int64,
|
|
|
|
n_bundles::Int64,
|
|
|
|
dir::String,
|
|
|
|
save_data::Bool,
|
2021-11-16 21:26:08 +00:00
|
|
|
)
|
2021-11-26 05:55:27 +00:00
|
|
|
bundle_snapshot_counter = 0
|
2021-11-16 21:26:08 +00:00
|
|
|
|
2021-12-07 03:19:35 +00:00
|
|
|
task::Union{Task,Nothing} = nothing
|
|
|
|
|
|
|
|
cl = CellList(args.particles_c, args.box; parallel=false)
|
|
|
|
|
2021-11-10 14:41:04 +00:00
|
|
|
start_time = now()
|
|
|
|
println("Started simulation at $start_time.")
|
|
|
|
|
2021-11-26 05:55:27 +00:00
|
|
|
@showprogress 0.6 for (integration_step, t) in enumerate(T0:δt:T)
|
|
|
|
if (integration_step % n_steps_before_snapshot == 0) && save_data
|
2021-11-26 13:50:47 +00:00
|
|
|
wait(task)
|
2021-11-26 05:55:27 +00:00
|
|
|
bundle_snapshot_counter += 1
|
|
|
|
save_snapshot!(args.bundle, bundle_snapshot_counter, t, args.particles)
|
|
|
|
|
|
|
|
if bundle_snapshot_counter == args.n_bundle_snapshots
|
2021-11-26 13:50:47 +00:00
|
|
|
task = @async begin
|
|
|
|
bundle_snapshot_counter = 0
|
|
|
|
n_bundles += 1
|
2021-11-26 05:55:27 +00:00
|
|
|
|
2021-11-26 13:50:47 +00:00
|
|
|
save_bundle(dir, args.bundle, n_bundles, t)
|
|
|
|
end
|
2021-11-26 05:55:27 +00:00
|
|
|
end
|
2021-11-10 14:41:04 +00:00
|
|
|
end
|
|
|
|
|
2021-11-10 23:24:07 +00:00
|
|
|
if integration_step % n_steps_before_verlet_list_update == 0
|
2021-12-07 03:19:35 +00:00
|
|
|
cl = update_verlet_list!(args, cl)
|
2021-11-10 14:41:04 +00:00
|
|
|
end
|
|
|
|
|
|
|
|
euler!(args)
|
|
|
|
end
|
|
|
|
|
2021-11-26 05:55:27 +00:00
|
|
|
if bundle_snapshot_counter > 0
|
|
|
|
bundle = first_n_snapshots(args.bundle, bundle_snapshot_counter)
|
|
|
|
|
|
|
|
n_bundles += 1
|
|
|
|
|
|
|
|
save_bundle(dir, bundle, n_bundles, T)
|
|
|
|
end
|
|
|
|
|
2021-11-10 14:41:04 +00:00
|
|
|
end_time = now()
|
2021-12-02 22:13:54 +00:00
|
|
|
elapsed_time = Dates.canonicalize(Dates.CompoundPeriod(end_time - start_time))
|
2021-11-18 18:42:11 +00:00
|
|
|
println("Simulation done at $end_time and took $elapsed_time.")
|
2021-11-10 14:41:04 +00:00
|
|
|
|
2021-11-26 05:55:27 +00:00
|
|
|
return nothing
|
2021-11-10 14:41:04 +00:00
|
|
|
end
|