diff --git a/Project.toml b/Project.toml index dd434a3..1c2f7c2 100644 --- a/Project.toml +++ b/Project.toml @@ -5,9 +5,11 @@ version = "0.2.0" [deps] BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" +CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b" CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0" CellListMap = "69e1c6dd-3888-40e6-b3c8-31ac5f578864" ColorSchemes = "35d6a980-a343-548e-a6ea-1d62b119f2f4" +DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" Dates = "ade2ca70-3891-5945-98fb-dc099432e06a" Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" Flux = "587475ba-b771-5e3f-ad9e-33799f191a9c" diff --git a/analysis/mean_squared_displacement.jl b/analysis/mean_squared_displacement.jl index 08348e8..98bec6a 100644 --- a/analysis/mean_squared_displacement.jl +++ b/analysis/mean_squared_displacement.jl @@ -9,6 +9,8 @@ using ProgressMeter: ProgressMeter using ReCo: ReCo +includet("../graphics/common_CairoMakie.jl") + # IMPORTANT: Disable the periodic boundary conditions # The arguments types have to match for the function to be overwritten! ReCo.push_to_verlet_list!(::Any, ::Any, ::Any) = nothing @@ -26,20 +28,6 @@ ReCo.minimum_image(v::SVector{2,Float64}, ::Float64) = v const δt = 1e-4 const Dₜ = ReCo.DEFAULT_Dₜ -function fill_with_bundle_property!( - v::Vector, property::Symbol, sim_dir::String, first_bundle::Int64=1 -) - bundle_paths = ReCo.sorted_bundle_paths(sim_dir) - - for i in first_bundle:length(bundle_paths) - bundle::ReCo.Bundle = JLD2.load_object(bundle_paths[i]) - - append!(v, getproperty(bundle, property)) - end - - return nothing -end - function max_possible_displacement(T::Float64, v₀::Float64, δt::Float64=δt, Dₜ::Float64=Dₜ) return T * v₀ + T / δt * sqrt(2 * Dₜ * δt) end @@ -52,6 +40,8 @@ function msd_simulation( parent_dir::String, comment::String="", ) + Random.seed!(42) + dir = ReCo.init_sim(; n_particles=1, v₀=v₀, @@ -72,10 +62,8 @@ function msd_simulation( end function mean_squared_displacement(; - n_simulations::Int64, v₀s::AbstractVector{Float64}, T::Float64 -) - Random.seed!(42) - + n_simulations::Int64, v₀s::NTuple{N,Float64}, T::Float64 +) where {N} n_v₀s = length(v₀s) main_parent_dir = "mean_squared_displacement_$(Dates.now())" @@ -99,7 +87,9 @@ function mean_squared_displacement(; end ts = Float64[] - fill_with_bundle_property!(ts, :t, sim_dirs[1, 1], 2) # Skip the first bundle to avoid t = 0 + ReCo.append_bundle_properties!( + (ts,), (:t,), sim_dirs[1, 1]; particle_slice=1, snapshot_slice=:, first_bundle=2 + ) # Skip the first bundle to avoid t = 0 mean_sq_displacements = zeros((length(ts), n_v₀s)) @@ -137,42 +127,9 @@ function expected_mean_squared_displacement(t::Float64, v₀::Float64) return (4 * Dₜ + 2 * v₀^2 / Dᵣ) * t + 2 * v₀^2 * (exp(-Dᵣ * t) - 1) / (Dᵣ^2) end -function init_cairomakie!() - CairoMakie.activate!() - set_theme!() - - return nothing -end - -function gen_figure() - text_width_in_pt = 405 - - return Figure(; - resolution=(text_width_in_pt, 0.55 * text_width_in_pt), - fontsize=10, - figure_padding=1, - ) -end - -function set_gaps!(fig::Figure) - colgap!(fig.layout, 5) - rowgap!(fig.layout, 5) - - return nothing -end - -function save_fig(filename::String, fig::Figure) - parent_dir = "exports/graphics" - mkpath(parent_dir) - - return save("$parent_dir/$filename", fig; pt_per_unit=1) -end - function plot_mean_sq_displacement_with_expectation( - ts::Vector{Float64}, - mean_sq_displacements::Matrix{Float64}, - v₀s::AbstractVector{Float64}, -) + ts::Vector{Float64}, mean_sq_displacements::Matrix{Float64}, v₀s::NTuple{N,Float64} +) where {N} init_cairomakie!() fig = gen_figure() @@ -205,7 +162,7 @@ function plot_mean_sq_displacement_with_expectation( end function run_msd_analysis() - v₀s = SVector(0.0, 20.0, 40.0, 60.0, 80.0) + v₀s = (0.0, 20.0, 40.0, 60.0, 80.0) ts, mean_sq_displacements = mean_squared_displacement(; n_simulations=200 * Threads.nthreads(), v₀s=v₀s, T=100.0 @@ -216,17 +173,17 @@ function run_msd_analysis() return nothing end -function plot_random_walk(T::Float64, v₀::Float64, seed::Int64) +function plot_random_walk(; T::Float64, v₀::Float64, seed::Int64) Random.seed!(seed) half_box_len = max_possible_displacement(T, v₀) dir = msd_simulation(v₀, half_box_len, T, 8.0, "random_walk_$(Dates.now())") ts = Float64[] - fill_with_bundle_property!(ts, :t, dir) - cs = SVector{2,Float64}[] - fill_with_bundle_property!(cs, :c, dir) + ReCo.append_bundle_properties!( + (ts, cs), (:t, :c), dir; particle_slice=1, snapshot_slice=: + ) init_cairomakie!() @@ -281,7 +238,7 @@ function plot_random_walk(T::Float64, v₀::Float64, seed::Int64) end function run_random_walk() - plot_random_walk(100_000.0, 0.0, 12345) + plot_random_walk(; T=100_000.0, v₀=0.0, seed=12345) return nothing end \ No newline at end of file diff --git a/analysis/pair_correlation.jl b/analysis/pair_correlation.jl deleted file mode 100644 index b19f952..0000000 --- a/analysis/pair_correlation.jl +++ /dev/null @@ -1,28 +0,0 @@ -if splitdir(pwd())[2] == "analysis" - cd("..") -end - -if splitdir(pwd())[2] != "ReCo.jl" - error("You have to be in the main directeory ReCo.jl!") -else - include("src/analysis/pair_correlation_function.jl") -end - -## - -using JLD2: JLD2 - -using CairoMakie -CairoMakie.activate!() -set_theme!(theme_black()) - -## - -dirs = readdir("exports/2021_11_19"; join=true) - -## - -for dir in dirs - data = JLD2.load("$dir/data.jld2") - display(pair_correlation(data["sol"], data["variables"]).fig) -end diff --git a/analysis/radial_distribution_function/g_of_r_d_ratio_with_0_45_packing_ratio.csv b/analysis/radial_distribution_function/g_of_r_d_ratio_with_0_45_packing_ratio.csv new file mode 100644 index 0000000..39928d0 --- /dev/null +++ b/analysis/radial_distribution_function/g_of_r_d_ratio_with_0_45_packing_ratio.csv @@ -0,0 +1,149 @@ +# Source: https://www.researchgate.net/publication/223333736_Reconstruction_of_Ross'_linear-mixing_model_for_the_equation_of_state_of_deuterium +r_d_ratio,g_of_r_d_ratio +0.99392,0.08604 +0.99388,0.18421 +0.99384,0.28847 +0.99382,0.34463 +0.99378,0.42994 +0.99372,0.52947 +0.99365,0.63441 +0.99361,0.68891 +0.99354,0.77406 +0.99344,0.87468 +0.99332,0.98013 +0.99325,1.03317 +0.99314,1.11822 +0.99297,1.21964 +0.99278,1.32539 +0.99267,1.37722 +0.99249,1.4622 +0.99224,1.56412 +0.99194,1.67 +0.99179,1.72089 +0.99152,1.8058 +0.99116,1.90793 +0.99074,2.01372 +0.99053,2.06399 +0.99015,2.1488 +0.98965,2.25086 +0.98909,2.35636 +0.98881,2.40637 +0.98831,2.49108 +0.98767,2.5925 +0.98699,2.69664 +0.98665,2.74669 +0.98607,2.83091 +0.98538,2.92885 +0.98467,3.02834 +0.98395,3.12761 +0.98362,3.17414 +0.98304,3.25677 +0.98242,3.34772 +0.98183,3.43776 +0.98129,3.52581 +0.98081,3.61118 +0.98043,3.69358 +0.98015,3.77322 +0.98,3.85156 +0.98014,3.95169 +0.98089,4.03823 +0.98357,4.10894 +0.9944,4.032 +1.00014,3.93792 +1.00441,3.85868 +1.00862,3.77571 +1.01272,3.69218 +1.01669,3.61072 +1.0216,3.51176 +1.02638,3.42062 +1.03128,3.33517 +1.03664,3.25085 +1.04281,3.16285 +1.04804,3.09295 +1.05614,2.98982 +1.06293,2.90588 +1.07039,2.81547 +1.07549,2.75441 +1.08352,2.65945 +1.09202,2.56047 +1.09758,2.49642 +1.10625,2.39796 +1.11214,2.33189 +1.12071,2.23742 +1.12939,2.14379 +1.13793,2.05435 +1.14612,1.97171 +1.15388,1.89756 +1.16359,1.81176 +1.17456,1.72585 +1.18615,1.64833 +1.19996,1.56268 +1.21169,1.47967 +1.22345,1.39143 +1.23543,1.30583 +1.24773,1.23043 +1.26375,1.14911 +1.28235,1.06205 +1.30026,0.97628 +1.32001,0.89731 +1.34607,0.8185 +1.38029,0.74139 +1.42373,0.67706 +1.47928,0.62965 +1.54985,0.64005 +1.60887,0.68208 +1.65956,0.742 +1.70693,0.80174 +1.75549,0.86677 +1.7997,0.93281 +1.8408,0.99879 +1.88278,1.07151 +1.92418,1.14479 +1.96069,1.20961 +2.00089,1.27998 +2.06624,1.29427 +2.12188,1.24296 +2.16454,1.1732 +2.20365,1.10578 +2.2468,1.03273 +2.28968,0.97094 +2.34711,0.92061 +2.40659,0.87706 +2.47617,0.86598 +2.54506,0.87775 +2.6091,0.91495 +2.67198,0.9548 +2.73415,0.99707 +2.79524,1.03776 +2.86004,1.06782 +2.92949,1.08716 +2.99726,1.09694 +3.06572,1.07401 +3.13305,1.0512 +3.20406,1.01975 +3.26722,0.99225 +3.33163,0.96833 +3.40415,0.95312 +3.46976,0.95561 +3.53944,0.96948 +3.61034,0.98695 +3.68191,1.00826 +3.74431,1.02656 +3.81351,1.04126 +3.88832,1.04634 +3.95948,1.04227 +4.02721,1.03357 +4.09307,1.02294 +4.16757,1.01022 +4.23256,0.99784 +4.30367,0.98757 +4.3786,0.98803 +4.44527,0.99118 +4.51289,0.99726 +4.58223,1.00774 +4.65288,1.01261 +4.72314,1.01386 +4.79174,1.02311 +4.86131,1.03134 +4.93191,1.0253 +4.98585,1.01943 diff --git a/analysis/radial_distribution_function/radial_distribution_function.jl b/analysis/radial_distribution_function/radial_distribution_function.jl new file mode 100644 index 0000000..5624d94 --- /dev/null +++ b/analysis/radial_distribution_function/radial_distribution_function.jl @@ -0,0 +1,216 @@ +using CairoMakie +using LaTeXStrings: @L_str +using StaticArrays: SVector +using JLD2: JLD2 +using Dates: Dates +using CSV: CSV +using DataFrames: DataFrames +using Random: Random + +using ReCo: ReCo + +includet("../../graphics/common_CairoMakie.jl") + +function radial_distribution_simulation(; + n_particles::Int64, v₀s::NTuple{N,Float64}, T::Float64, packing_ratio::Float64 +) where {N} + Random.seed!(42) + + n_v₀s = length(v₀s) + + sim_dirs = Vector{String}(undef, n_v₀s) + + parent_dir = "radial_distribution_$(Dates.now())" + + Threads.@threads for v₀_ind in 1:n_v₀s + v₀ = v₀s[v₀_ind] + + dir = ReCo.init_sim(; + n_particles=n_particles, + v₀=v₀, + packing_ratio=packing_ratio, + parent_dir=parent_dir, + comment="$v₀", + ) + + ReCo.run_sim(dir; duration=T, seed=v₀_ind) + + sim_dirs[v₀_ind] = dir + end + + return sim_dirs +end + +function circular_shell_volume(lower_radius, Δradius) + return π * (2 * lower_radius * Δradius + Δradius^2) +end + +function radial_distribution(; + sim_dirs::Vector{String}, n_radii::Int64, n_last_snapshots::Int64, n_particles::Int64 +) + sim_consts::ReCo.SimConsts = JLD2.load_object("$(sim_dirs[1])/sim_consts.jld2") + + particle_radius = sim_consts.particle_radius + + min_lower_radius = 0.0 + max_lower_radius = 5 * (2 * particle_radius) + Δradius = (max_lower_radius - min_lower_radius) / n_radii + + lower_radii = LinRange(min_lower_radius, max_lower_radius, n_radii) + + box_volume = (2 * sim_consts.half_box_len)^2 + volume_per_particle = box_volume / n_particles + + n_sim_dirs = length(sim_dirs) + + gs = Vector{Vector{Float64}}(undef, n_sim_dirs) + + Threads.@threads for sim_dir_ind in 1:n_sim_dirs + sim_dir = sim_dirs[sim_dir_ind] + + cs = Matrix{SVector{2,Float64}}(undef, (n_particles, n_last_snapshots)) + + bundle_paths = ReCo.sorted_bundle_paths(sim_dir; rev=true) + + snapshot_conunter = 0 + break_bundle_path_loop = false + + for bundle_path in bundle_paths + bundle::ReCo.Bundle = JLD2.load_object(bundle_path) + + for snapshot_ind in (bundle.n_snapshots):-1:1 + snapshot_conunter += 1 + + @simd for particle_ind in 1:n_particles + cs[particle_ind, snapshot_conunter] = bundle.c[ + particle_ind, snapshot_ind + ] + end + + if snapshot_conunter == n_last_snapshots + break_bundle_path_loop = true + break + end + end + + if break_bundle_path_loop + break + end + end + + if snapshot_conunter != n_last_snapshots + error("snapshot_conunter != n_last_snapshots") + end + + g = zeros(n_radii) + + for snapshot_ind in 1:n_last_snapshots + for p1_ind in 1:n_particles + for p2_ind in (p1_ind + 1):n_particles + c1 = cs[p1_ind, snapshot_ind] + c2 = cs[p2_ind, snapshot_ind] + + r⃗₁₂ = c2 - c1 + r⃗₁₂ = ReCo.minimum_image(r⃗₁₂, sim_consts.half_box_len) + + distance = ReCo.norm2d(r⃗₁₂) + + lower_radius_ind = ceil(Int64, distance / Δradius) + + if lower_radius_ind <= n_radii + g[lower_radius_ind] += 2 + end + end + end + end + + for (lower_radius_ind, lower_radius) in enumerate(lower_radii) + g[lower_radius_ind] *= + volume_per_particle / ( + n_last_snapshots * + n_particles * + circular_shell_volume(lower_radius, Δradius) + ) + end + + gs[sim_dir_ind] = g + end + + return (lower_radii, gs, particle_radius) +end + +function plot_radial_distributions( + v₀s::NTuple{N,Float64}, + lower_radii::AbstractVector{Float64}, + gs::Vector{Vector{Float64}}, + particle_radius::Float64, +) where {N} + println("Plotting the radial distributions") + + init_cairomakie!() + + fig = gen_figure() + + max_lower_radius = maximum(lower_radii) + max_g = maximum(maximum.(gs)) + + ax = Axis( + fig[1:2, 1:2]; + xticks=0:(2 * particle_radius):ceil(Int64, max_lower_radius), + yticks=0:ceil(Int64, max_g), + xlabel=L"r / d", + ylabel=L"g", + ) + + lines!( + ax, + SVector(0.0, max_lower_radius), + SVector(1.0, 1.0); + linestyle=:dash, + color=:red, + linewidth=1, + ) + + expected = CSV.read( + "analysis/g_of_r_d_ratio_with_0_45_packing_ratio.csv", + DataFrames.DataFrame; + header=2, + ) + lines!( + ax, + expected.r_d_ratio, + expected.g_of_r_d_ratio; + label="Expectation for dense hard spheres", + color=:green, + ) + + for (g_ind, g) in enumerate(gs) + scatterlines!( + ax, lower_radii, g; markersize=3, linewidth=0.5, label=L"v_0 = %$(v₀s[g_ind])" + ) + end + + axislegend(ax; position=:rt, padding=3, rowgap=-3) + + save_fig("radial_distribution.pdf", fig) + + return nothing +end + +function run_radial_distribution_analysis() + v₀s = (0.0, 80.0) + + n_particles = 1000 + + sim_dirs = radial_distribution_simulation(; + n_particles=n_particles, v₀s=v₀s, T=100.0, packing_ratio=0.45 + ) + + lower_radii, gs, particle_radius = radial_distribution(; + sim_dirs=sim_dirs, n_radii=75, n_last_snapshots=200, n_particles=n_particles + ) + + plot_radial_distributions(v₀s, lower_radii, gs, particle_radius) + + return nothing +end \ No newline at end of file diff --git a/graphics/common_CairoMakie.jl b/graphics/common_CairoMakie.jl new file mode 100644 index 0000000..3862278 --- /dev/null +++ b/graphics/common_CairoMakie.jl @@ -0,0 +1,29 @@ +function init_cairomakie!() + CairoMakie.activate!() + set_theme!() + + return nothing +end + +function gen_figure(; padding=4) + text_width_in_pt = 405 + + return Figure(; + resolution=(text_width_in_pt, 0.55 * text_width_in_pt), + fontsize=10, + figure_padding=padding, + ) +end + +function set_gaps!(fig::Figure) + colgap!(fig.layout, 5) + rowgap!(fig.layout, 5) + + return nothing +end + +function save_fig(filename::String, fig::Figure, parent_dir="exports/graphics") + mkpath(parent_dir) + + return save("$parent_dir/$filename", fig; pt_per_unit=1) +end \ No newline at end of file diff --git a/graphics/pontential.jl b/graphics/pontential.jl index 3de23b4..60c2fd2 100644 --- a/graphics/pontential.jl +++ b/graphics/pontential.jl @@ -1,6 +1,8 @@ using CairoMakie using LaTeXStrings: @L_str +includet("common_CairoMakie.jl") + const minimum_r_σ_ratio = 2^(1 / 6) function U_LJ_ϵ_ratio(r_σ_ratio::Real) @@ -16,36 +18,41 @@ function U_WCA_ϵ_ratio(r_σ_ratio::Real) end end -text_width_in_pt = 405 +function plot_potentials() + init_cairomakie!() -fig = Figure(; - resolution=(text_width_in_pt, 0.55 * text_width_in_pt), fontsize=10, figure_padding=1 -) + fig = gen_figure() -max_x = 2.5 + max_x = 2.5 -max_y = 1.05 -min_y = -max_y + max_y = 1.05 + min_y = -max_y -ax = Axis(fig[1, 1]; xlabel=L"r/σ", ylabel=L"U/ϵ", limits=(0.88, max_x, min_y, max_y)) + ax = Axis(fig[1, 1]; xlabel=L"r/σ", ylabel=L"U/ϵ", limits=(0.88, max_x, min_y, max_y)) -r_σ_ratio = LinRange(0.95, max_x, 1000) + r_σ_ratio = LinRange(0.95, max_x, 1000) -LJ = lines!(ax, r_σ_ratio, U_LJ_ϵ_ratio.(r_σ_ratio)) + LJ = lines!(ax, r_σ_ratio, U_LJ_ϵ_ratio.(r_σ_ratio)) -WCA = lines!(ax, r_σ_ratio, U_WCA_ϵ_ratio.(r_σ_ratio)) + WCA = lines!(ax, r_σ_ratio, U_WCA_ϵ_ratio.(r_σ_ratio)) -minimum_r_σ_ratio_line = lines!( - ax, [minimum_r_σ_ratio, minimum_r_σ_ratio], [min_y, max_y]; linestyle=:dash, linewidth=1 -) + minimum_r_σ_ratio_line = lines!( + ax, + [minimum_r_σ_ratio, minimum_r_σ_ratio], + [min_y, max_y]; + linestyle=:dash, + linewidth=1, + ) -Legend( - fig[1, 2], - [LJ, WCA, minimum_r_σ_ratio_line], - [L"U_{LJ}", L"U_{WCA}", L"\frac{r}{σ} = 2^{1/6}"], -) + Legend( + fig[1, 2], + [LJ, WCA, minimum_r_σ_ratio_line], + [L"U_{LJ}", L"U_{WCA}", L"\frac{r}{σ} = 2^{1/6}"], + ) -colgap!(fig.layout, 5) -rowgap!(fig.layout, 5) + set_gaps!(fig) -save("exports/graphics/potential.pdf", fig; pt_per_unit=1) \ No newline at end of file + save_fig("potential.pdf", fig) + + return nothing +end \ No newline at end of file diff --git a/src/RL/RL.jl b/src/RL/RL.jl index 419acee..857d079 100644 --- a/src/RL/RL.jl +++ b/src/RL/RL.jl @@ -124,12 +124,7 @@ function run_rl(; agent(PRE_EPISODE_STAGE, env) # Episode - ReCo.run_sim( - dir; - duration=episode_duration, - seed=rand(1:typemax(Int64)), - env_helper=env_helper, - ) + ReCo.run_sim(dir; duration=episode_duration, seed=episode, env_helper=env_helper) env.shared.terminated = true diff --git a/src/analysis/pair_correlation_function.jl b/src/analysis/pair_correlation_function.jl deleted file mode 100644 index 7d23ed1..0000000 --- a/src/analysis/pair_correlation_function.jl +++ /dev/null @@ -1,83 +0,0 @@ -using CairoMakie, LaTeXStrings -using LoopVectorization: @turbo - -using ReCo: minimum_image, norm2d - -function plot_g(radius, g, variables) - fig = Figure() - ax = Axis( - fig[1, 1]; - xticks=0:ceil(Int64, maximum(radius)), - yticks=0:ceil(Int64, maximum(g)), - xlabel=L"r", - ylabel=L"g(r)", - title="v₀ = $(variables.v₀)", - ) - - scatterlines!(ax, radius, g; color=:white, markercolor=:red) - - return fig -end - -function pair_correlation(sol, variables) - n_r = 100 - n_last_frames = 200 - min_radius = variables.particle_diameter / 2 - max_radius = 10.0 * variables.particle_diameter - dr = (max_radius - min_radius) / n_r - - radius = range(min_radius, max_radius; length=n_r) - N_g = zeros(variables.n_particles, n_r) - - Threads.@threads for r_ind in 1:n_r - r = radius[r_ind] - - @simd for i in 1:(variables.n_particles) - for j in 1:(variables.n_particles) - if i != j - for k in 1:n_last_frames - frame = variables.n_snapshots - k + 1 - c1 = sol.center[i, frame] - c2 = sol.center[j, frame] - - r⃗₁₂ = c1 - c2 - - r⃗₁₂ = minimum_image(r⃗₁₂, variables.half_box_len) - - distance = norm2d(r⃗₁₂) - - if (distance >= r) && (distance <= r + dr) - N_g[i, r_ind] += 1 - end - end - end - end - end - end - - g = zeros(n_r) - - @simd for r_ind in 1:n_r - r = radius[r_ind] - tmp_g = 0.0 - - @turbo for i in 1:(variables.n_particles) - tmp_g += N_g[i, r_ind] - end - - tmp_g *= - (2 * variables.half_box_len)^2 / ( - (variables.n_particles * n_last_frames) * - variables.n_particles * - 2 * - π * - r * - dr - ) - g[r_ind] = tmp_g - end - - fig = plot_g(radius, g, variables) - - return (; fig, radius, g) -end diff --git a/src/data.jl b/src/data.jl index 7e2e284..ff13118 100644 --- a/src/data.jl +++ b/src/data.jl @@ -39,6 +39,9 @@ struct Bundle t::Vector{Float64} c::Matrix{SVector{2,Float64}} φ::Matrix{Float64} + + n_particles::Int64 + n_snapshots::Int64 end function Bundle(n_particles::Int64, n_snapshots::Int64) @@ -46,11 +49,13 @@ function Bundle(n_particles::Int64, n_snapshots::Int64) Vector{Float64}(undef, n_snapshots), Matrix{SVector{2,Float64}}(undef, (n_particles, n_snapshots)), Matrix{Float64}(undef, (n_particles, n_snapshots)), + n_particles, + n_snapshots, ) end function first_n_snapshots(bundle::Bundle, n::Int64) - return Bundle(bundle.t[1:n], bundle.c[:, 1:n], bundle.φ[:, 1:n]) + return Bundle(bundle.t[1:n], bundle.c[:, 1:n], bundle.φ[:, 1:n], bundle.n_particles, n) end function save_snapshot!( @@ -78,7 +83,7 @@ function save_bundle(dir::String, bundle::Bundle, n_bundle::Int64, T::Float64) return nothing end -function sorted_bundle_paths(dir::String) +function sorted_bundle_paths(dir::String; rev::Bool=false) bundle_paths = readdir("$dir/bundles"; join=true, sort=false) n_bundles = length(bundle_paths) @@ -93,7 +98,37 @@ function sorted_bundle_paths(dir::String) bundle_nums[i] = parse(Int64, bundle_num_string) end - sort_perm = sortperm(bundle_nums) + sort_perm = sortperm(bundle_nums; rev=rev) return bundle_paths[sort_perm] +end + +function append_bundle_properties!( + vs::NTuple{N,Vector}, + properties::NTuple{N,Symbol}, + sim_dir::String; + particle_slice=nothing, + snapshot_slice, + first_bundle::Int64=1, +) where {N} + bundle_paths = ReCo.sorted_bundle_paths(sim_dir) + + for i in first_bundle:length(bundle_paths) + bundle::ReCo.Bundle = JLD2.load_object(bundle_paths[i]) + + for (v_ind, v) in enumerat(vs) + property = properties[v_ind] + bundle_property = getproperty(bundle, property) + + if !isnothing(particle_slice) + bundle_property_view = view(bundle_property, particle_slice, snapshot_slice) + else + bundle_property_view = view(bundle_property, snapshot_slice) + end + + append!(v, bundle_property_view) + end + end + + return nothing end \ No newline at end of file