const DEFAULT_PACKING_RATIO = 0.5 const DEFAULT_δt = 1e-5 const DEFAULT_SKIN_TO_INTERACTION_R_RATIO = 2.0 const DEFAULT_EXPORTS_DIR = "exports" const DEFAULT_PARENT_DIR = "" const DEFAULT_COMMENT = "" const DEFAULT_PARTICLE_RADIUS = 0.5 const DEFAULT_Dₜ = 1.0 const DEFAULT_μₜ = 1.0 const DEFAULT_ϵ = 100.0 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 gen_particles(grid_n::Int64, grid_box_width::Float64, half_box_len::Float64) particles = Vector{Particle}(undef, grid_n^2) id = 1 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, half_box_len), rand(Uniform(-π, π)), ) id += 1 end end return particles end function gen_particles(bundle::Bundle, n_particles::Int64) particles = Vector{Particle}(undef, n_particles) @simd for id in 1:n_particles particles[id] = Particle(id, bundle.c[id, end], bundle.φ[id, end]) end return particles end function gen_sim_consts( n_particles::Int64, v₀::Float64; δt::Float64=DEFAULT_δt, packing_ratio::Float64=DEFAULT_PACKING_RATIO, skin_to_interaction_r_ratio::Float64=DEFAULT_SKIN_TO_INTERACTION_R_RATIO, half_box_len::Float64=0.0, ) @assert n_particles > 0 @assert v₀ >= 0 @assert δt in 1e-7:1e-7:1e-4 @assert packing_ratio > 0 μₜ = DEFAULT_μₜ Dₜ = DEFAULT_Dₜ particle_radius = DEFAULT_PARTICLE_RADIUS Dᵣ = 3 * Dₜ / ((2 * particle_radius)^2) σ = 2 * particle_radius * 2^(-1 / 6) ϵ = DEFAULT_ϵ interaction_r = 2^(1 / 6) * σ skin_r = skin_to_interaction_r_ratio * interaction_r buffer = 2.5 max_approach_after_one_integration_step = buffer * (2 * v₀ * δt) @assert skin_r >= interaction_r + max_approach_after_one_integration_step if v₀ != 0.0 n_steps_before_verlet_list_update = floor( Int64, (skin_r - interaction_r) / max_approach_after_one_integration_step ) else n_steps_before_verlet_list_update = 1000 end grid_n = round(Int64, ceil(sqrt(n_particles))) n_particles = grid_n^2 if half_box_len == 0.0 half_box_len = sqrt(n_particles * π / packing_ratio) * σ / 4 elseif packing_ratio != DEFAULT_PACKING_RATIO error("You can not specify half_box_len and packing_ratio at the same time!") end grid_box_width = 2 * half_box_len / grid_n return SimConsts( # Input n_particles, v₀, δt, packing_ratio, # Internal μₜ, Dₜ, particle_radius, Dᵣ, σ, ϵ, interaction_r, skin_r, n_steps_before_verlet_list_update, grid_n, half_box_len, grid_box_width, ) end function init_sim_with_sim_consts( sim_consts::SimConsts; exports_dir::String=DEFAULT_EXPORTS_DIR, parent_dir::String=DEFAULT_PARENT_DIR, comment::String=DEFAULT_COMMENT, ) particles = gen_particles( sim_consts.grid_n, sim_consts.grid_box_width, sim_consts.half_box_len ) bundle = Bundle(sim_consts.n_particles, 1) save_snapshot!(bundle, 1, 0.0, particles) dir = exports_dir if length(parent_dir) > 0 dir *= "/$parent_dir" end start_datetime = Dates.now() dir *= "/$(start_datetime)_N=$(sim_consts.n_particles)_v=$(sim_consts.v₀)_#$(rand(1000:9999))" if length(comment) > 0 dir *= "_$comment" end mkpath(dir) task = @async JLD2.save_object("$dir/sim_consts.jld2", sim_consts) save_bundle(dir, bundle, 1, 0.0) runs_dir = "$dir/runs" mkpath(runs_dir) wait(task) return dir end function init_sim(; n_particles::Int64, v₀::Float64, δt::Float64=DEFAULT_δt, packing_ratio::Float64=DEFAULT_PACKING_RATIO, skin_to_interaction_r_ratio::Float64=DEFAULT_SKIN_TO_INTERACTION_R_RATIO, exports_dir::String=DEFAULT_EXPORTS_DIR, parent_dir::String=DEFAULT_PARENT_DIR, comment::String=DEFAULT_COMMENT, half_box_len::Float64=0.0, ) sim_consts = gen_sim_consts( n_particles, v₀; δt=δt, packing_ratio=packing_ratio, skin_to_interaction_r_ratio=skin_to_interaction_r_ratio, half_box_len, ) return init_sim_with_sim_consts( sim_consts; exports_dir=exports_dir, parent_dir=parent_dir, comment=comment ) end