From 95d3d5d9afe67845116be5bd4c360a3852a31dba Mon Sep 17 00:00:00 2001 From: MoBit Date: Sat, 20 Nov 2021 14:15:22 +0100 Subject: [PATCH] pair_correlation --- .gitignore | 4 +- Project.toml | 1 + analysis/pair_correlation.jl | 26 ++++++++ src/ReCo.jl | 4 +- src/analysis/pair_correlation_function.jl | 72 +++++++++++++++++++++++ src/animation.jl | 38 +++++++----- src/pair_correlation_function.jl | 67 --------------------- 7 files changed, 129 insertions(+), 83 deletions(-) create mode 100644 analysis/pair_correlation.jl create mode 100644 src/analysis/pair_correlation_function.jl delete mode 100644 src/pair_correlation_function.jl diff --git a/.gitignore b/.gitignore index ccd7092..11ad1ca 100644 --- a/.gitignore +++ b/.gitignore @@ -1,2 +1,4 @@ Manifest.toml -exports \ No newline at end of file +exports +.ipynb_checkpoints +*.ipynb \ No newline at end of file diff --git a/Project.toml b/Project.toml index 53a08f7..684253a 100644 --- a/Project.toml +++ b/Project.toml @@ -1,5 +1,6 @@ [deps] BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" +CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0" ColorSchemes = "35d6a980-a343-548e-a6ea-1d62b119f2f4" Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" GLMakie = "e9467ef8-e4e7-5192-8a1a-b1aee30e663a" diff --git a/analysis/pair_correlation.jl b/analysis/pair_correlation.jl new file mode 100644 index 0000000..a78c43c --- /dev/null +++ b/analysis/pair_correlation.jl @@ -0,0 +1,26 @@ +if splitdir(pwd())[2] == "analysis" + cd("..") +end + +if splitdir(pwd())[2] != "ReCo" + error("You have to be in the main directeory ReCo!") +else + include("src/ReCo.jl") + includet("src/analysis/pair_correlation_function.jl") +end + +## + +using CairoMakie +CairoMakie.activate!() + +## + +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/src/ReCo.jl b/src/ReCo.jl index e8927bb..d630980 100644 --- a/src/ReCo.jl +++ b/src/ReCo.jl @@ -16,6 +16,8 @@ import Base: push!, run, iterate using Revise using BenchmarkTools +set_theme!(theme_black()) + includet("PreVector.jl") includet("Particle.jl") includet("setup.jl") @@ -23,4 +25,4 @@ includet("simulation.jl") includet("data.jl") includet("animation.jl") includet("postprocess.jl") -includet("run.jl") +includet("run.jl") \ No newline at end of file diff --git a/src/analysis/pair_correlation_function.jl b/src/analysis/pair_correlation_function.jl new file mode 100644 index 0000000..95d6b8e --- /dev/null +++ b/src/analysis/pair_correlation_function.jl @@ -0,0 +1,72 @@ +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, n_r) + + Threads.@threads for r_ind in 1:n_r + r = radius[r_ind] + + @simd for i in 1:(variables.N) + for j in 1:(variables.N) + if i != j + for k in 1:n_last_frames + frame = variables.n_frames - k + 1 + c1 = sol.center[i, frame] + c2 = sol.center[j, frame] + + r⃗₁₂ = SVector(c1[1] - c2[1], c1[2] - c2[2]) + + r⃗₁₂ = minimum_image(r⃗₁₂; l=variables.l) + + distance = sqrt(r⃗₁₂[1]^2 + r⃗₁₂[2]^2) + + 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 + + @tturbo for i in 1:(variables.N) + tmp_g += N_g[i, r_ind] + end + + tmp_g *= + (2 * variables.l)^2 / + ((variables.N * n_last_frames) * variables.N * 2 * π * r * dr) + g[r_ind] = tmp_g + end + + fig = plot_g(radius, g, variables) + + return (; fig, radius, g) +end diff --git a/src/animation.jl b/src/animation.jl index 3bfde1f..d463cbf 100644 --- a/src/animation.jl +++ b/src/animation.jl @@ -1,8 +1,6 @@ -function animate(directory::String, sol::Solution, variables) +function animate(directory::String, sol::Solution, variables; framerate=0) println("Generating animation...") - set_theme!(theme_black()) - fig = Figure(; resolution=(1080, 1080)) ax = Axis( fig[1, 1]; @@ -18,7 +16,14 @@ function animate(directory::String, sol::Solution, variables) animation_path = "$directory/animation.mkv" - record(fig, animation_path; framerate=variables.framerate) do io + record_framerate = 0 + if framerate > 0 + record_framerate = framerate + else + record_framerate = variables.framerate + end + + record(fig, animation_path; framerate=record_framerate) do io circles = Observable(Vector{Circle}(undef, variables.N)) colors = Observable(Vector{RGBAf}(undef, variables.N)) @@ -33,13 +38,14 @@ function animate(directory::String, sol::Solution, variables) skin_colors = Observable(Vector{RGBAf}(undef, variables.N)) end + π2 = 2 * π + particle_radius = variables.particle_diameter / 2 + @showprogress 0.6 for frame in 1:(variables.n_frames) @simd for i in 1:(variables.N) - circles[][i] = Circle( - Point2(sol.center[i, frame]), variables.particle_diameter / 2 - ) + circles[][i] = Circle(Point2(sol.center[i, frame]), particle_radius) - color = get(color_scheme, rem2pi(sol.φ[i, frame] / (2 * π), RoundDown)) + color = get(color_scheme, rem2pi(sol.φ[i, frame], RoundDown) / π2) colors[][i] = RGBAf(color) if variables.debug @@ -50,19 +56,23 @@ function animate(directory::String, sol::Solution, variables) Point2(sol.center[i, frame]), variables.skin_r ) - interaction_colors[][i] = RGBAf(color, 0.12) - skin_colors[][i] = RGBAf(color, 0.06) + interaction_colors[][i] = RGBAf(color, 0.08) + skin_colors[][i] = RGBAf(color, 0.04) end end if frame > 1 if variables.debug @simd for i in 1:(variables.N) - segments_x[][2 * i - 1] = sol.center[i, frame - 1][1] - segments_x[][2 * i] = sol.center[i, frame][1] + first_ind = 2 * i - 1 + second_ind = 2 * i + frame_min_1 = frame - 1 - segments_y[][2 * i - 1] = sol.center[i, frame - 1][2] - segments_y[][2 * i] = sol.center[i, frame][2] + segments_x[][first_ind] = sol.center[i, frame_min_1][1] + segments_x[][second_ind] = sol.center[i, frame][1] + + segments_y[][first_ind] = sol.center[i, frame_min_1][2] + segments_y[][second_ind] = sol.center[i, frame][2] end if frame == 2 diff --git a/src/pair_correlation_function.jl b/src/pair_correlation_function.jl deleted file mode 100644 index 4bcd1ec..0000000 --- a/src/pair_correlation_function.jl +++ /dev/null @@ -1,67 +0,0 @@ -function pair_correlation(frame=0, dr=0.1) - sol, args, filename = run(; - N=1000, - T=1.0, - v=20.0, - δt=1e-5, - save_at=0.1, - framerate=10, - n_steps_before_verlet_list_update=1000, - ) - - n_r = 100 - radius = range(0, args.l; length=n_r) - N_g = zeros((args.N, n_r)) - - if frame == 0 - frame = args.n_frames - end - - for r_ind in 1:n_r - r = radius[r_ind] - - for i in 1:(args.N) - for j in 1:(args.N) - if i != j - c1 = sol.center[i, frame] - c2 = sol.center[j, frame] - - r⃗₁₂ = SVector{2}(c1) - SVector{2}(c2) - - r⃗₁₂ = minimum_image(r⃗₁₂; l=args.l) - - distance = sqrt(r⃗₁₂[1]^2 + r⃗₁₂[2]^2) - - if (distance >= r) && (distance <= r + dr) - N_g[i, r_ind] += 1 - end - end - end - end - end - - g = zeros(n_r) - - for r_ind in 1:n_r - r = radius[r_ind] - tmp_g = 0.0 - - for i in 1:(args.N) - tmp_g += N_g[i, r_ind] - end - - tmp_g = tmp_g * (2 * args.l)^2 / (args.N^2 * 2 * π * r * dr) - g[r_ind] = tmp_g - end - - return plot_g(radius, g) -end - -function plot_g(radius, g) - fig = Figure(; resolution=(1080, 1080)) - ax = Axis(fig[1, 1]; xlabel=L"r", ylabel=L"g(r)") - - scatterlines!(ax, radius, g) - - return (; fig, radius, g) -end