From f4734064c5e359f90455db9e3a5c79a307f46c2d Mon Sep 17 00:00:00 2001 From: MoBit Date: Tue, 16 Nov 2021 22:26:08 +0100 Subject: [PATCH] pair_correlation_funciton --- src/PreVector.jl | 8 ++++ src/ReCo.jl | 2 +- src/animation.jl | 13 ++----- src/pair_correlation_function.jl | 67 ++++++++++++++++++++++++++++++++ src/run.jl | 45 ++++++++++----------- src/simulation.jl | 26 +++++++------ 6 files changed, 118 insertions(+), 43 deletions(-) create mode 100644 src/pair_correlation_function.jl diff --git a/src/PreVector.jl b/src/PreVector.jl index 20f2919..a0a94ee 100644 --- a/src/PreVector.jl +++ b/src/PreVector.jl @@ -17,3 +17,11 @@ function reset!(pv::PreVector{T}) where {T} return nothing end + +function iterate(pv::PreVector{T}, state=1) where {T} + if state > pv.last_ind + return nothing + else + return (pv.v[state], state + 1) + end +end \ No newline at end of file diff --git a/src/ReCo.jl b/src/ReCo.jl index 1db082f..d83f0ff 100644 --- a/src/ReCo.jl +++ b/src/ReCo.jl @@ -8,7 +8,7 @@ using GLMakie, ColorSchemes, LaTeXStrings using StaticArrays import Dates: now, CompoundPeriod, canonicalize -import Base: push!, run +import Base: push!, run, iterate # Development deps using Revise diff --git a/src/animation.jl b/src/animation.jl index a11b64b..a1cadbb 100644 --- a/src/animation.jl +++ b/src/animation.jl @@ -12,12 +12,9 @@ function animate(sol::Solution, args, name_part::String; framerate::Int64=10) ylabel=L"y", ) - Colorbar( - fig[1, 2]; - limits=(0, 2), - colormap=ColorSchemes.cyclic_mrybm_35_75_c68_n256_s25, - label=L"\frac{\varphi}{\pi}", - ) + color_scheme = ColorSchemes.cyclic_mrybm_35_75_c68_n256_s25 + + Colorbar(fig[1, 2]; limits=(0, 2), colormap=color_scheme, label=L"\frac{\varphi}{\pi}") animation_path = "exports/$name_part.mkv" @@ -42,9 +39,7 @@ function animate(sol::Solution, args, name_part::String; framerate::Int64=10) Point2(sol.center[i, frame]), args.particle_diameter / 2 ) - color = get( - ColorSchemes.rainbow, rem2pi(sol.φ[i, frame] / (2 * π), RoundDown) - ) + color = get(color_scheme, rem2pi(sol.φ[i, frame] / (2 * π), RoundDown)) colors[][i] = RGBAf(color) if args.debug diff --git a/src/pair_correlation_function.jl b/src/pair_correlation_function.jl new file mode 100644 index 0000000..fbb6196 --- /dev/null +++ b/src/pair_correlation_function.jl @@ -0,0 +1,67 @@ +function pair_correlation(frame=0, dr=0.1) + sol, args, name_part = 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 + rs = 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 = rs[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 = rs[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(rs, g) +end + +function plot_g(rs, g) + fig = Figure(; resolution=(1080, 1080)) + ax = Axis(fig[1, 1]; xlabel=L"r", ylabel=L"g(r)") + + scatterlines!(ax, rs, g) + + return (fig, rs, g) +end \ No newline at end of file diff --git a/src/run.jl b/src/run.jl index 0ce476e..e6ed9be 100644 --- a/src/run.jl +++ b/src/run.jl @@ -8,8 +8,7 @@ function run(; save_data::Bool=false, n_steps_before_verlet_list_update::Int64=100, debug::Bool=false, - ) - +) Random.seed!(42) μ = 1.0 @@ -33,27 +32,29 @@ function run(; n_steps_before_save = round(Int64, save_at / δt) 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), - skin_r = skin_r, - skin_r² = skin_r^2, - verlet_list = [PreVector(Int64, N - 1) for i in 1:(N - 1)], - n_frames = floor(Int64, integration_steps / n_steps_before_save) + 1, - debug = debug, + 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), + skin_r=skin_r, + skin_r²=skin_r^2, + verlet_list=[PreVector(Int64, N - 1) for i in 1:(N - 1)], + n_frames=floor(Int64, integration_steps / n_steps_before_save) + 1, + debug=debug, ) - sol, end_time = simulate(args, δt, T, n_steps_before_verlet_list_update, n_steps_before_save) + sol, end_time = simulate( + args, δt, T, n_steps_before_verlet_list_update, n_steps_before_save + ) name_part = "$(end_time)_T=$(T)_N=$(N)_v=$(v)_dt=$(δt)" @@ -65,5 +66,5 @@ function run(; animate(sol, args, name_part; framerate=framerate) end - return (sol = sol, args = args, name_part = name_part) + return (sol=sol, args=args, name_part=name_part) end diff --git a/src/simulation.jl b/src/simulation.jl index cc74a18..be0599f 100644 --- a/src/simulation.jl +++ b/src/simulation.jl @@ -4,9 +4,9 @@ function update_verlet_list!(args) @simd for pv in args.verlet_list reset!(pv) end - + for i in 1:(args.N - 1) - for j in (i + 1):args.N + for j in (i + 1):(args.N) p1 = args.particles[i] p2 = args.particles[j] @@ -17,6 +17,8 @@ function update_verlet_list!(args) end end end + + return nothing end function euler!(args) @@ -24,17 +26,19 @@ function euler!(args) p = args.particles[i] verlet_list = args.verlet_list[p.id] - for j in 1:verlet_list.last_ind + for j in 1:(verlet_list.last_ind) p2 = args.particles[verlet_list.v[j]] - overlapping, r⃗₁₂, distance² = are_overlapping(p, p2, args.interaction_r², args.l) + overlapping, r⃗₁₂, distance² = are_overlapping( + p, p2, args.interaction_r², args.l + ) if overlapping c = args.c₁ / (distance²^4) * (args.c₂ / (distance²^3) - 1) - + @simd for k in 1:2 dck = c * r⃗₁₂[k] - + p.tmp_c[k] -= dck p2.tmp_c[k] += dck end @@ -59,19 +63,19 @@ function euler!(args) return nothing end -function simulate(args, +function simulate( + args, δt::Float64, T::Float64, n_steps_before_verlet_list_update::Int64, n_steps_before_save::Int64, - ) - +) sol = Solution(args.N, args.n_frames) frame = 0 - + frame = save_frame!(sol, frame, 0.0, args.particles) - + start_time = now() println("Started simulation at $start_time.")