1
0
Fork 0
mirror of https://gitlab.rlp.net/mobitar/ReCo.jl.git synced 2024-12-29 17:53:31 +00:00

Bundles for resuming simulations

This commit is contained in:
MoBit 2021-11-26 06:55:27 +01:00
parent 10018160e2
commit e2f06265df
8 changed files with 208 additions and 151 deletions

View file

@ -24,6 +24,5 @@ includet("Particle.jl")
includet("setup.jl")
includet("simulation.jl")
includet("data.jl")
includet("animation.jl")
includet("postprocess.jl")
# includet("animation.jl")
includet("run.jl")

View file

@ -31,7 +31,7 @@ function pair_correlation(sol, variables)
for j in 1:(variables.N)
if i != j
for k in 1:n_last_frames
frame = variables.n_frames - k + 1
frame = variables.n_snapshots - k + 1
c1 = sol.center[i, frame]
c2 = sol.center[j, frame]

View file

@ -41,7 +41,7 @@ function animate(directory::String, sol::Solution, variables; framerate=0)
π2 = 2 * π
particle_radius = variables.particle_diameter / 2
@showprogress 0.6 for frame in 1:(variables.n_frames)
@showprogress 0.6 for frame in 1:(variables.n_snapshots)
@simd for i in 1:(variables.N)
circles[][i] = Circle(Point2(sol.center[i, frame]), particle_radius)

View file

@ -1,52 +1,60 @@
struct Solution
struct Bundle
t::Vector{Float64}
center::Matrix{Vector{Float64}}
c::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_frame!(sol, frame, t, particles)
frame += 1
function Bundle(N::Int64, n_snapshots::Int64)
return Bundle(
zeros(n_snapshots),
[zeros(2) for i in 1:N, j in 1:n_snapshots],
zeros(N, n_snapshots),
)
end
sol.t[frame] = t
function first_n_snapshots(bundle::Bundle, n::Int64)
return Bundle(bundle.t[1:n], bundle.c[:, 1:n], bundle.φ[:, 1:n])
end
function save_snapshot!(
bundle::Bundle, n_snapshot::Int64, t::Float64, particles::Vector{Particle}
)
bundle.t[n_snapshot] = t
bundle_c = bundle.c
bundle_φ = bundle.φ
@simd for p in particles
p_id = p.id
p_c = p.c
pre_pos = sol.center[p_id, frame]
bundle_c_id = bundle_c[p_id, n_snapshot]
@turbo for i in 1:2
pre_pos[i] = p.c[i]
@turbo for j in 1:2
bundle_c_id[j] = p_c[j]
end
sol.φ[p_id, frame] = p.φ
bundle_φ[p_id, n_snapshot] = p.φ
end
return frame
end
function save_data_jld(directory::String; kwargs...)
path = "$directory/data.jld2"
JLD2.jldsave(path; kwargs...)
println("Data saved to $path.")
return nothing
end
function save_variables(directory::String; kwargs...)
path = "$directory/variables.json"
open(path, "w") do io
JSON3.pretty(io, kwargs...)
function set_status(dir::String, n_bundles::Int64, T::Float64)
open("$dir/status.json", "w") do f
JSON3.write(f, Dict("n_bundles" => n_bundles, "T" => round(T; digits=3)))
end
println("Variables saved to $path.")
return nothing
end
function save_bundle(dir::String, bundle::Bundle, n::Int64, T::Float64)
bundle_dir = "$dir/bundles/bundle$n"
mkpath(bundle_dir)
JLD2.save_object("$bundle_dir/bundle.jld2", bundle)
set_status(dir, n, T)
return nothing
end

View file

@ -1,23 +0,0 @@
function postprocess(sol, variables, save_data)
if !isdir("exports")
mkdir("exports")
end
directory = "exports/$(variables.end_time)_T=$(variables.T)_N=$(variables.N)_v=$(variables.v)"
if save_data || (variables.framerate > 0)
mkdir(directory)
save_variables(directory; variables=variables)
end
if save_data
save_data_jld(directory; sol=sol, variables=variables)
end
if variables.framerate > 0
animate(directory, sol, variables)
end
return nothing
end

View file

@ -1,108 +1,110 @@
function run(;
N::Int64,
T::Int64,
v::Float64,
δt::Float64=1e-5,
function run_sim(
dir::String;
duration::Float64,
snapshot_at::Float64=0.1,
framerate::Int64=0,
save_data::Bool=true,
n_steps_before_verlet_list_update::Int64=1000,
packing_ratio::Float64=0.5,
debug::Bool=false,
seed::Int64=42,
n_bundle_snapshots::Int64=100,
)
@assert N > 0
@assert T > 0
@assert v >= 0
@assert δt > 0
@assert snapshot_at 0.01:0.01:T
@assert framerate >= 0
@assert length(dir) > 0
@assert duration > 0
@assert snapshot_at in 0.01:0.01:duration
@assert n_steps_before_verlet_list_update >= 0
@assert packing_ratio > 0
@assert seed > 0
@assert n_bundle_snapshots >= 0
Random.seed!(seed)
μ = 1.0
D₀ = 1.0
particle_diameter = 1.0
Dᵣ = 3 * D₀ / (particle_diameter^2)
σ = 1.0
ϵ = 100.0
interaction_r = 2^(1 / 6) * σ
sim_consts = JSON3.read(read("$dir/sim_consts.json", String))
grid_n = round(Int64, ceil(sqrt(N)))
skin_r =
1.1 * (
2 * sim_consts.v * n_steps_before_verlet_list_update * sim_consts.δt +
2 * sim_consts.interaction_r
)
N = grid_n^2
integration_steps = floor(Int64, duration / sim_consts.δt) + 1
n_steps_before_snapshot = round(Int64, snapshot_at / sim_consts.δt)
l = sqrt(N * π / packing_ratio) * σ / 4
grid_box_width = 2 * l / grid_n
n_snapshots = floor(Int64, integration_steps / n_steps_before_snapshot) + 1
skin_r = 1.1 * (2 * v * n_steps_before_verlet_list_update * δt + 2 * interaction_r)
integration_steps = floor(Int64, T / δt) + 1
n_steps_before_save = round(Int64, snapshot_at / δt)
n_frames = floor(Int64, integration_steps / n_steps_before_save) + 1
n_bundle_snapshots = min(n_snapshots, n_bundle_snapshots)
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),
v=sim_consts.v,
skin_r=skin_r,
skin_r²=skin_r^2,
verlet_list=[PreVector(Int64, N - 1) for i in 1:(N - 1)],
n_frames=n_frames,
debug=debug,
n_snapshots=n_snapshots,
c₁=4 * sim_consts.ϵ * 6 * sim_consts.σ^6 * sim_consts.δt * sim_consts.μ,
c₂=2 * sim_consts.σ^6,
c₃=sqrt(2 * sim_consts.D₀ * sim_consts.δt),
c₄=sqrt(2 * sim_consts.Dᵣ * sim_consts.δt),
vδt=sim_consts.v * sim_consts.δt,
μ=sim_consts.μ,
interaction_r=sim_consts.interaction_r,
interaction_r²=sim_consts.interaction_r^2,
N=sim_consts.N,
l=sim_consts.l,
particle_diameter=sim_consts.particle_diameter,
particles=generate_particles(
sim_consts.grid_n, sim_consts.grid_box_width, sim_consts.l
),
verlet_list=[PreVector{Int64}(sim_consts.N - 1) for i in 1:(sim_consts.N - 1)],
n_bundle_snapshots=n_bundle_snapshots,
bundle=Bundle(sim_consts.N, n_bundle_snapshots),
)
sol, end_time = simulate(
args, δt, T, n_steps_before_verlet_list_update, n_steps_before_save
)
status = JSON3.read(read("$dir/status.json", String))
T0::Float64 = status.T
args = nothing
T = T0 + duration
variables = (;
start_datetime = now()
run_vars = (;
# Input
N,
T,
v,
δt,
duration,
snapshot_at,
framerate,
n_steps_before_verlet_list_update,
packing_ratio,
debug,
seed,
n_bundle_snapshots,
# Calculated
μ,
D₀,
particle_diameter,
Dᵣ,
σ,
ϵ,
interaction_r,
grid_n,
l,
grid_box_width,
skin_r,
integration_steps,
n_steps_before_save,
n_frames,
# Output
end_time,
n_steps_before_snapshot,
n_snapshots,
T,
# Read
T0,
start_datetime,
)
postprocess(sol, variables, save_data)
n_bundles::Int64 = status.n_bundles
next_bundle = n_bundles + 1
next_bundle_path = "$dir/bundles/bundle$next_bundle"
mkpath(next_bundle_path)
return (; sol, variables)
if n_bundle_snapshots > 0
save_data = true
open("$next_bundle_path/run_vars.json", "w") do f
JSON3.pretty(f, run_vars)
end
else
save_data = false
end
simulate(
args,
sim_consts.δt,
T0,
T,
n_steps_before_verlet_list_update,
n_steps_before_snapshot,
n_bundles,
dir,
save_data,
)
return dir
end

View file

@ -1,14 +1,11 @@
function initial_particle_grid_pos(i, j; grid_box_width, l)
dim1_pos(x) = (x - 0.5) * grid_box_width - l
return SVector(dim1_pos(i), dim1_pos(j))
term = -0.5 * grid_box_width - l
return [k * grid_box_width + term for k in (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
@ -25,3 +22,61 @@ function generate_particles(grid_n, grid_box_width, l)
return particles
end
function init_sim(;
N::Int64, v::Float64, δt::Float64=1e-5, packing_ratio::Float64=0.5, comment::String=""
)
@assert N > 0
@assert v >= 0
@assert δt in 1e-7:1e-7:1e-4
@assert packing_ratio > 0
μ = 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 = sqrt(N * π / packing_ratio) * σ / 4
grid_box_width = 2 * l / grid_n
sim_consts = (;
# Input
N,
v,
δt,
packing_ratio,
# Calculated
μ,
D₀,
particle_diameter,
Dᵣ,
σ,
ϵ,
interaction_r,
grid_n,
l,
grid_box_width,
)
particles = generate_particles(grid_n, grid_box_width, l)
bundle = Bundle(N, 1)
save_snapshot!(bundle, 1, 0.0, particles)
dir = "exports/$(now())_N=$(N)_v=$(v)_#$(rand(1000:9999))$comment"
mkpath(dir)
open("$dir/sim_consts.json", "w") do f
JSON3.pretty(f, sim_consts)
end
save_bundle(dir, bundle, 1, 0.0)
return dir
end

View file

@ -66,22 +66,30 @@ end
function simulate(
args,
δt::Float64,
T::Int64,
T0::Float64,
T::Float64,
n_steps_before_verlet_list_update::Int64,
n_steps_before_save::Int64,
n_steps_before_snapshot::Int64,
n_bundles::Int64,
dir::String,
save_data::Bool,
)
sol = Solution(args.N, args.n_frames)
frame = 0
frame = save_frame!(sol, frame, 0.0, args.particles)
bundle_snapshot_counter = 0
start_time = now()
println("Started simulation at $start_time.")
@showprogress 0.6 for (integration_step, t) in enumerate(0:δt:T)
if integration_step % n_steps_before_save == 0
frame = save_frame!(sol, frame, t, args.particles)
@showprogress 0.6 for (integration_step, t) in enumerate(T0:δt:T)
if (integration_step % n_steps_before_snapshot == 0) && save_data
bundle_snapshot_counter += 1
save_snapshot!(args.bundle, bundle_snapshot_counter, t, args.particles)
if bundle_snapshot_counter == args.n_bundle_snapshots
bundle_snapshot_counter = 0
n_bundles += 1
save_bundle(dir, args.bundle, n_bundles, t)
end
end
if integration_step % n_steps_before_verlet_list_update == 0
@ -91,9 +99,17 @@ function simulate(
euler!(args)
end
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
end_time = now()
elapsed_time = canonicalize(CompoundPeriod(end_time - start_time))
println("Simulation done at $end_time and took $elapsed_time.")
return (; sol, end_time)
return nothing
end