2021-12-02 22:13:54 +00:00
|
|
|
using ProgressMeter: @showprogress
|
|
|
|
using LoopVectorization: @turbo
|
|
|
|
using Distributions: Normal
|
|
|
|
using Dates: Dates, now
|
|
|
|
import Base: wait
|
|
|
|
|
2021-11-10 14:41:04 +00:00
|
|
|
rand_normal01() = rand(Normal(0, 1))
|
|
|
|
|
|
|
|
function update_verlet_list!(args)
|
|
|
|
@simd for pv in args.verlet_list
|
|
|
|
reset!(pv)
|
|
|
|
end
|
2021-11-16 21:26:08 +00:00
|
|
|
|
2021-11-10 14:41:04 +00:00
|
|
|
for i in 1:(args.N - 1)
|
2021-11-16 21:26:08 +00:00
|
|
|
for j in (i + 1):(args.N)
|
2021-11-10 14:41:04 +00:00
|
|
|
p1 = args.particles[i]
|
|
|
|
p2 = args.particles[j]
|
|
|
|
|
2021-11-17 18:47:05 +00:00
|
|
|
overlapping = are_overlapping(p1, p2, args.skin_r², args.l).overlapping
|
2021-11-10 14:41:04 +00:00
|
|
|
|
|
|
|
if overlapping
|
|
|
|
push!(args.verlet_list[i], j)
|
|
|
|
end
|
|
|
|
end
|
|
|
|
end
|
2021-11-16 21:26:08 +00:00
|
|
|
|
|
|
|
return nothing
|
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)
|
|
|
|
p = args.particles[i]
|
2021-11-10 14:41:04 +00:00
|
|
|
verlet_list = args.verlet_list[p.id]
|
|
|
|
|
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(
|
|
|
|
p, p2, args.interaction_r², args.l
|
|
|
|
)
|
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-11-26 02:35:39 +00:00
|
|
|
p.tmp_c[k] -= dc
|
|
|
|
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()
|
|
|
|
|
|
|
|
restrict_coordinates!(p; l=args.l)
|
|
|
|
|
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-11-10 14:41:04 +00:00
|
|
|
start_time = now()
|
|
|
|
println("Started simulation at $start_time.")
|
|
|
|
|
2021-11-26 13:50:47 +00:00
|
|
|
task::Union{Task,Nothing} = nothing
|
|
|
|
|
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-11-10 14:41:04 +00:00
|
|
|
update_verlet_list!(args)
|
|
|
|
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
|