From 389f99221b855f29465113b4437179f7b343e308 Mon Sep 17 00:00:00 2001 From: MoBit Date: Tue, 7 Dec 2021 05:52:10 +0100 Subject: [PATCH] SVector --- src/Particle.jl | 25 +++++++---------------- src/analysis/pair_correlation_function.jl | 6 +++--- src/animation.jl | 3 ++- src/data.jl | 23 ++++++--------------- src/run.jl | 7 +------ src/setup.jl | 4 ++-- src/shape.jl | 9 +++----- src/simulation.jl | 25 +++++++++-------------- 8 files changed, 34 insertions(+), 68 deletions(-) diff --git a/src/Particle.jl b/src/Particle.jl index 272fce0..83e46fd 100644 --- a/src/Particle.jl +++ b/src/Particle.jl @@ -1,16 +1,15 @@ using StaticArrays: SVector -using LoopVectorization: @turbo mutable struct Particle id::Int64 - c::Vector{Float64} # Center - tmp_c::Vector{Float64} # Temporary center + c::SVector{2,Float64} # Center + tmp_c::SVector{2,Float64} # Temporary center φ::Float64 # Angle - function Particle(id::Int64, c::Vector{Float64}, φ::Float64) - return new(id, c, copy(c), φ) + function Particle(id::Int64, c::SVector{2,Float64}, φ::Float64) + return new(id, c, c, φ) end end @@ -25,9 +24,7 @@ function restrict_coordinate(value::Float64, l::Float64) end function restrict_coordinates!(p::Particle, l::Float64) - @simd for i in 1:2 - p.tmp_c[i] = restrict_coordinate(p.tmp_c[i], l) - end + p.tmp_c = SVector{2}(restrict_coordinate(p.tmp_c[i], l) for i in 1:2) return nothing end @@ -43,11 +40,11 @@ function minimum_image_coordinate(value::Float64, l::Float64) end function minimum_image(v::SVector{2,Float64}, l::Float64) - return SVector(minimum_image_coordinate(v[1], l), minimum_image_coordinate(v[2], l)) + return SVector{2}(minimum_image_coordinate(v[i], l) for i in 1:2) end function are_overlapping(p1::Particle, p2::Particle, overlapping_r²::Float64, l::Float64) - r⃗₁₂ = SVector(p2.c[1] - p1.c[1], p2.c[2] - p1.c[2]) # 1 -> 2 + r⃗₁₂ = p2.c - p1.c # 1 -> 2 r⃗₁₂ = minimum_image(r⃗₁₂, l) @@ -57,11 +54,3 @@ function are_overlapping(p1::Particle, p2::Particle, overlapping_r²::Float64, l return (; overlapping, r⃗₁₂, distance²) end - -function update!(p::Particle) - @turbo for i in 1:2 - p.c[i] = p.tmp_c[i] - end - - return nothing -end diff --git a/src/analysis/pair_correlation_function.jl b/src/analysis/pair_correlation_function.jl index 6a6f842..7b0765e 100644 --- a/src/analysis/pair_correlation_function.jl +++ b/src/analysis/pair_correlation_function.jl @@ -1,6 +1,6 @@ using CairoMakie, LaTeXStrings using StaticArrays -using LoopVectorization: @tturbo +using LoopVectorization: @turbo using ReCo: minimum_image @@ -41,7 +41,7 @@ function pair_correlation(sol, variables) c1 = sol.center[i, frame] c2 = sol.center[j, frame] - r⃗₁₂ = SVector(c1[1] - c2[1], c1[2] - c2[2]) + r⃗₁₂ = c1 - c2 r⃗₁₂ = minimum_image(r⃗₁₂, variables.l) @@ -62,7 +62,7 @@ function pair_correlation(sol, variables) r = radius[r_ind] tmp_g = 0.0 - @tturbo for i in 1:(variables.N) + @turbo for i in 1:(variables.N) tmp_g += N_g[i, r_ind] end diff --git a/src/animation.jl b/src/animation.jl index 1657591..e6b6aa5 100644 --- a/src/animation.jl +++ b/src/animation.jl @@ -20,7 +20,8 @@ function animate_bundle!(args) for frame in 1:length(bundle_t) @simd for i in 1:(args.N) - args.circles[][i] = Circle(Point2(bundle_c[i, frame]), args.particle_radius) + c = bundle_c[i, frame] + args.circles[][i] = Circle(Point2(c[1], c[2]), args.particle_radius) color = get(args.color_scheme, rem2pi(bundle_φ[i, frame], RoundDown) / π2) args.colors[][i] = RGBAf(color) diff --git a/src/data.jl b/src/data.jl index ecd2bf6..514bb21 100644 --- a/src/data.jl +++ b/src/data.jl @@ -1,18 +1,17 @@ using StaticArrays: SVector -using LoopVectorization: @turbo using JLD2: JLD2 struct Bundle t::Vector{Float64} - c::Matrix{Vector{Float64}} + c::Matrix{SVector{2,Float64}} φ::Matrix{Float64} end 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), + Vector{Float64}(undef, n_snapshots), + Matrix{SVector{2,Float64}}(undef, (N, n_snapshots)), + Matrix{Float64}(undef, (N, n_snapshots)), ) end @@ -25,20 +24,10 @@ function save_snapshot!( ) bundle.t[n_snapshot] = t - bundle_c = bundle.c - bundle_φ = bundle.φ - @simd for p in particles - p_id = p.id - p_c = p.c + bundle.c[p.id, n_snapshot] = p.c - bundle_c_id = bundle_c[p_id, n_snapshot] - - @turbo for j in 1:2 - bundle_c_id[j] = p_c[j] - end - - bundle_φ[p_id, n_snapshot] = p.φ + bundle.φ[p.id, n_snapshot] = p.φ end return nothing diff --git a/src/run.jl b/src/run.jl index 570b245..f67b87e 100644 --- a/src/run.jl +++ b/src/run.jl @@ -43,11 +43,6 @@ function run_sim( particles = generate_particles(bundle, sim_consts.N) - particles_c = [SVector(0.0, 0.0) for i in 1:(sim_consts.N)] - for i in 1:(sim_consts.N) - particles_c[i] = SVector(particles[i].c[1], particles[i].c[2]) - end - args = ( v=sim_consts.v, skin_r=skin_r, @@ -65,7 +60,7 @@ function run_sim( l=sim_consts.l, particle_diameter=sim_consts.particle_diameter, particles=particles, - particles_c=particles_c, + particles_c=[particles[i].c for i in 1:(sim_consts.N)], 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), diff --git a/src/setup.jl b/src/setup.jl index 2e6684d..ef27e91 100644 --- a/src/setup.jl +++ b/src/setup.jl @@ -4,7 +4,7 @@ using JSON3: JSON3 function initial_particle_grid_pos(i::Int64, j::Int64, grid_box_width::Float64, l::Float64) term = -0.5 * grid_box_width - l - return [k * grid_box_width + term for k in (i, j)] + return SVector{2}(k * grid_box_width + term for k in (i, j)) end function generate_particles(grid_n::Int64, grid_box_width::Float64, l::Float64) @@ -28,7 +28,7 @@ end function generate_particles(bundle::Bundle, N::Int64) particles = Vector{Particle}(undef, N) - for id in 1:N + @simd for id in 1:N particles[id] = Particle(id, bundle.c[id, end], bundle.φ[id, end]) end diff --git a/src/shape.jl b/src/shape.jl index 4dc816f..8bcc67b 100644 --- a/src/shape.jl +++ b/src/shape.jl @@ -12,8 +12,8 @@ function project_back_from_unit_circle(θ, l) end function center_of_mass(args) - x_proj_sum = SVector(0, 0) - y_proj_sum = SVector(0, 0) + x_proj_sum = SVector(0.0, 0.0) + y_proj_sum = SVector(0.0, 0.0) for p in args.particles x_proj_sum += project_to_unit_circle(p.c[1], args.l) @@ -37,10 +37,7 @@ function center_of_mass(args) y_θ = atan(y_proj_sum[2], y_proj_sum[1]) end - return SVector( - project_back_from_unit_circle(x_θ, args.l), - project_back_from_unit_circle(y_θ, args.l), - ) + return SVector{2}(project_back_from_unit_circle(θ, args.l) for θ in (x_θ, y_θ)) end function gyration_tensor(args) diff --git a/src/simulation.jl b/src/simulation.jl index e301c3a..c66074a 100644 --- a/src/simulation.jl +++ b/src/simulation.jl @@ -1,8 +1,8 @@ using Dates: Dates, now using ProgressMeter: @showprogress -using LoopVectorization: @turbo using Distributions: Normal using CellListMap: Box, CellList, map_pairwise!, UpdateCellList! +using StaticArrays: SVector import Base: wait @@ -23,8 +23,8 @@ function update_verlet_list!(args, cl) reset!(pv) end - for i in 1:(args.N) - args.particles_c[i] = SVector(args.particles[i].c[1], args.particles[i].c[2]) + @simd for i in 1:(args.N) + args.particles_c[i] = args.particles[i].c end cl = UpdateCellList!(args.particles_c, args.box, cl; parallel=false) @@ -55,28 +55,23 @@ function euler!(args) if overlapping c = args.c₁ / (distance²^4) * (args.c₂ / (distance²^3) - 1) - @turbo for k in 1:2 - dc = c * r⃗₁₂[k] - - p1.tmp_c[k] -= dc - p2.tmp_c[k] += dc - end + dc = c * r⃗₁₂ + p1.tmp_c -= dc + p2.tmp_c += dc end end end @simd for p in args.particles - e = SVector(cos(p.φ), sin(p.φ)) - - @turbo for i in 1:2 - p.tmp_c[i] += args.vδt * e[i] + args.c₃ * rand_normal01() - end + p.tmp_c += + args.vδt * SVector(cos(p.φ), sin(p.φ)) + + args.c₃ * SVector(rand_normal01(), rand_normal01()) p.φ += args.c₄ * rand_normal01() restrict_coordinates!(p, args.l) - update!(p) + p.c = p.tmp_c end return nothing