diff --git a/.gitignore b/.gitignore index 11ad1ca..e2dc786 100644 --- a/.gitignore +++ b/.gitignore @@ -1,4 +1,9 @@ Manifest.toml exports + +# Jupyter notebooks .ipynb_checkpoints -*.ipynb \ No newline at end of file +*.ipynb + +# Profiling +*.mem \ No newline at end of file diff --git a/graphics/center_of_mass.jl b/graphics/center_of_mass.jl index fd37944..1e0d33c 100644 --- a/graphics/center_of_mass.jl +++ b/graphics/center_of_mass.jl @@ -29,7 +29,8 @@ M = R * ((Ap + Bp) / 2) θ = atan(-M[2], M[1]) COM = project_back_from_unit_circle(θ, L) -COMp = R * Point(cos(θ), -sin(θ)) +si, co = sincos(θ) +COMp = R * Point(co, -si) ## diff --git a/src/Particle.jl b/src/Particle.jl index 83e46fd..712b420 100644 --- a/src/Particle.jl +++ b/src/Particle.jl @@ -23,8 +23,12 @@ function restrict_coordinate(value::Float64, l::Float64) return value end +function restrict_coordinates!(v::SVector{2,Float64}, l::Float64) + return SVector(restrict_coordinate(v[1], l), restrict_coordinate(v[2], l)) +end + function restrict_coordinates!(p::Particle, l::Float64) - p.tmp_c = SVector{2}(restrict_coordinate(p.tmp_c[i], l) for i in 1:2) + p.tmp_c = restrict_coordinates!(p.tmp_c, l) return nothing end @@ -40,11 +44,13 @@ function minimum_image_coordinate(value::Float64, l::Float64) end function minimum_image(v::SVector{2,Float64}, l::Float64) - return SVector{2}(minimum_image_coordinate(v[i], l) for i in 1:2) + return SVector(minimum_image_coordinate(v[1], l), minimum_image_coordinate(v[2], l)) end -function are_overlapping(p1::Particle, p2::Particle, overlapping_r²::Float64, l::Float64) - r⃗₁₂ = p2.c - p1.c # 1 -> 2 +function are_overlapping( + c1::SVector{2,Float64}, c2::SVector{2,Float64}, overlapping_r²::Float64, l::Float64 +) + r⃗₁₂ = c2 - c1 # 1 -> 2 r⃗₁₂ = minimum_image(r⃗₁₂, l) diff --git a/src/setup.jl b/src/setup.jl index ef27e91..a0f4d65 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 SVector{2}(k * grid_box_width + term for k in (i, j)) + return SVector(i * grid_box_width + term, j * grid_box_width + term) end function generate_particles(grid_n::Int64, grid_box_width::Float64, l::Float64) diff --git a/src/shape.jl b/src/shape.jl index 8bcc67b..5019125 100644 --- a/src/shape.jl +++ b/src/shape.jl @@ -1,23 +1,24 @@ using StaticArrays: SVector, SMatrix using LinearAlgebra: eigvals, Hermitian -function project_to_unit_circle(x, l) +function project_to_unit_circle(x::Float64, l::Float64) φ = (x + l) * π / l - return SVector(cos(φ), sin(φ)) + si, co = sincos(φ) + return SVector(co, si) end -function project_back_from_unit_circle(θ, l) +function project_back_from_unit_circle(θ::T, l::Float64) where {T<:Real} x = θ * l / π - l return restrict_coordinate(x, l) end -function center_of_mass(args) +function center_of_mass(particles::Vector{Particle}, l::Float64) 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) - y_proj_sum += project_to_unit_circle(p.c[2], args.l) + for p in particles + x_proj_sum += project_to_unit_circle(p.c[1], l) + y_proj_sum += project_to_unit_circle(p.c[2], l) end # Prevent for example atan(1e-16, 1e-15) != 0 with rounding @@ -37,30 +38,31 @@ function center_of_mass(args) y_θ = atan(y_proj_sum[2], y_proj_sum[1]) end - return SVector{2}(project_back_from_unit_circle(θ, args.l) for θ in (x_θ, y_θ)) + COM_x = project_back_from_unit_circle(x_θ, l) + COM_y = project_back_from_unit_circle(y_θ, l) + + return SVector(COM_x, COM_y) end -function gyration_tensor(args) - COM = center_of_mass(args) +function gyration_tensor(particles::Vector{Particle}, l::Float64) + COM = center_of_mass(particles, l) S11 = 0.0 S12 = 0.0 S22 = 0.0 - for p in args.particles - c = p.c - x = restrict_coordinate(c[1] - COM[1], args.l) - y = restrict_coordinate(c[2] - COM[2], args.l) + for p in particles + shifted_c = restrict_coordinates!(p.c - COM, l) - S11 += x^2 - S12 += x * y - S22 += y^2 + S11 += shifted_c[1]^2 + S12 += shifted_c[1] * shifted_c[2] + S22 += shifted_c[2]^2 end return Hermitian(SMatrix{2,2}(S11, S12, S12, S22)) end -function gyration_tensor_eigvals_ratio(args) - ev = eigvals(gyration_tensor(args)) # Eigenvalues are sorted +function gyration_tensor_eigvals_ratio(particles::Vector{Particle}, l::Float64) + ev = eigvals(gyration_tensor(particles, l)) # Eigenvalues are sorted return ev[1] / ev[2] end \ No newline at end of file diff --git a/src/simulation.jl b/src/simulation.jl index 1584881..add6857 100644 --- a/src/simulation.jl +++ b/src/simulation.jl @@ -41,19 +41,20 @@ function update_verlet_list!(args, cl) end function euler!(args) - for i in 1:(args.N - 1) - p1 = args.particles[i] - verlet_list = args.verlet_list[p1.id] + for id1 in 1:(args.N - 1) + p1 = args.particles[id1] + p1_c = p1.c + verlet_list = args.verlet_list[id1] - for j in 1:(verlet_list.last_ind) - p2 = args.particles[verlet_list.v[j]] + for id2 in view(verlet_list.v, 1:(verlet_list.last_ind)) + p2 = args.particles[id2] overlapping, r⃗₁₂, distance² = are_overlapping( - p1, p2, args.interaction_r², args.l + p1_c, p2.c, args.interaction_r², args.l ) if overlapping - c = args.c₁ / (distance²^4) * (args.c₂ / (distance²^3) - 1) + c = args.c₁ / (distance²^4) * (args.c₂ / (distance²^3) - 1.0) dc = c * r⃗₁₂ p1.tmp_c -= dc @@ -63,9 +64,11 @@ function euler!(args) end @simd for p in args.particles - p.tmp_c += - args.vδt * SVector(cos(p.φ), sin(p.φ)) + - args.c₃ * SVector(rand_normal01(), rand_normal01()) + si, co = sincos(p.φ) + p.tmp_c += SVector( + args.vδt * co + args.c₃ * rand_normal01(), + args.vδt * si + args.c₃ * rand_normal01(), + ) p.φ += args.c₄ * rand_normal01()