Optimizations

.gitignore

@@ -1,4 +1,9 @@
 Manifest.toml
 exports
+
+# Jupyter notebooks
 .ipynb_checkpoints
 *.ipynb
+
+# Profiling
+*.mem

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)
 
 ##
 

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)
+function are_overlapping(
-    r⃗₁₂ = p2.c - p1.c # 1 -> 2
+    c1::SVector{2,Float64}, c2::SVector{2,Float64}, overlapping_r²::Float64, l::Float64
+)
+    r⃗₁₂ = c2 - c1 # 1 -> 2
 
     r⃗₁₂ = minimum_image(r⃗₁₂, l)
 

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)

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
+    for p in particles
-        x_proj_sum += project_to_unit_circle(p.c[1], args.l)
+        x_proj_sum += project_to_unit_circle(p.c[1], l)
-        y_proj_sum += project_to_unit_circle(p.c[2], args.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)
+function gyration_tensor(particles::Vector{Particle}, l::Float64)
-    COM = center_of_mass(args)
+    COM = center_of_mass(particles, l)
 
     S11 = 0.0
     S12 = 0.0
     S22 = 0.0
 
-    for p in args.particles
+    for p in particles
-        c = p.c
+        shifted_c = restrict_coordinates!(p.c - COM, l)
-        x = restrict_coordinate(c[1] - COM[1], args.l)
-        y = restrict_coordinate(c[2] - COM[2], args.l)
 
-        S11 += x^2
+        S11 += shifted_c[1]^2
-        S12 += x * y
+        S12 += shifted_c[1] * shifted_c[2]
-        S22 += y^2
+        S22 += shifted_c[2]^2
     end
 
     return Hermitian(SMatrix{2,2}(S11, S12, S12, S22))
 end
 
-function gyration_tensor_eigvals_ratio(args)
+function gyration_tensor_eigvals_ratio(particles::Vector{Particle}, l::Float64)
-    ev = eigvals(gyration_tensor(args)) # Eigenvalues are sorted
+    ev = eigvals(gyration_tensor(particles, l)) # Eigenvalues are sorted
     return ev[1] / ev[2]
 end

src/simulation.jl

@@ -41,19 +41,20 @@ function update_verlet_list!(args, cl)
 end
 
 function euler!(args)
-    for i in 1:(args.N - 1)
+    for id1 in 1:(args.N - 1)
-        p1 = args.particles[i]
+        p1 = args.particles[id1]
-        verlet_list = args.verlet_list[p1.id]
+        p1_c = p1.c
+        verlet_list = args.verlet_list[id1]
 
-        for j in 1:(verlet_list.last_ind)
+        for id2 in view(verlet_list.v, 1:(verlet_list.last_ind))
-            p2 = args.particles[verlet_list.v[j]]
+            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 +=
+        si, co = sincos(p.φ)
-            args.vδt * SVector(cos(p.φ), sin(p.φ)) +
+        p.tmp_c += SVector(
-            args.c₃ * SVector(rand_normal01(), rand_normal01())
+            args.vδt * co + args.c₃ * rand_normal01(),
+            args.vδt * si + args.c₃ * rand_normal01(),
+        )
 
         p.φ += args.c₄ * rand_normal01()