diff --git a/src/boxgrid.jl b/src/boxgrid.jl index 8118365394281ae870f3e98c6a447b2b03d38b6c..ec9cfa9e0d1dd0110916c8be66901b5d041bb708 100644 --- a/src/boxgrid.jl +++ b/src/boxgrid.jl @@ -158,3 +158,119 @@ function relevant_boxes(so::AbstractObstacle, bg::BoxGrid, w::AbstractDomain{3}; @warn "Automatic box detection for 3D obstacles is not supported. Override `relevant_boxes` for your obstacle to manually select boxes" maxlog=1 return vec(CartesianIndices(bg.boxes)) end + + + +## BoxGrid pairwise interaction computation + +@forcedef pairwise_interactions(interaction::Function, sim::Simulation, bg::BoxGrid) = + pairwise_interactions(interaction, sim, bg, Val{isthreaded(sim)}()) + + +""" + pairwise_interactions(interaction::Function, sim::Simulation, bg::BoxGrid, ::Val{false}) +Uses a [`BoxGrid`](@ref) to efficiently compute pairwise interactions. +`BoxGrid` structure is mutated in the sense that particles are sorted into their boxes. +""" +@forcedef function pairwise_interactions(interaction::Function, sim::Simulation, bg::BoxGrid, ::Val{false}) + # For execution order preservation compared to multithreaded simulations the + # boxes are iterated in the same striped manner. + boxparticles!(bg, particles(sim)) + + numstripes, numrows... = size(bg.boxes) + rowit = CartesianIndices(numrows) + + for i = 1:2:numstripes + for j = rowit + I = CartesianIndex(i,j) + isempty(bg.boxes[I]) || @forcefun pairwise_interactions(interaction, sim, bg.boxes[I]) + for J ∈ bg.neighborlist[I] + if !isempty(bg.boxes[J]) + @forcefun pairwise_interactions(interaction, sim, bg.boxes[I], bg.boxes[J]) + end + end # for J ∈ neighbors + end + end + for i = 2:2:numstripes + for j = rowit + I = CartesianIndex(i,j) + isempty(bg.boxes[I]) || @forcefun pairwise_interactions(interaction, sim, bg.boxes[I]) + for J ∈ bg.neighborlist[I] + if !isempty(bg.boxes[J]) + @forcefun pairwise_interactions(interaction, sim, bg.boxes[I], bg.boxes[J]) + end + end # for J ∈ neighbors + end + end +end + + +""" + pairwise_interactions(interaction::Function, sim::Simulation, bg::BoxGrid, ::Val{true}) +Multi-threaded version of [`pairwise_interactions!_serial`](@ref). + +## Extended help + +Multithreading is implemented by dividing the system into two sets of interleaved stripes (slabs in 3D) +of boxes, each stripe being one box wide. +This division ensures that each stripe is independent of all other stripes in same set. + +Therefore, we can first integrate all the stripes in the first set in parallel, synchronise, then +integrate the second set. +""" +@forcedef function pairwise_interactions(interaction::Function, sim::Simulation, bg::BoxGrid, ::Val{true}) + boxparticles!(bg, particles(sim)) + + numstripes, numrows... = size(bg.boxes) + rowit = CartesianIndices(numrows) + + Threads.@threads for i = 1:2:numstripes + for j = rowit + I = CartesianIndex(i,j) + isempty(bg.boxes[I]) || @forcefun pairwise_interactions(interaction, sim, bg.boxes[I]) + for J ∈ bg.neighborlist[I] + if !isempty(bg.boxes[J]) + @forcefun pairwise_interactions(interaction, sim, bg.boxes[I], bg.boxes[J]) + end + end # for J ∈ neighbors + end + end + Threads.@threads for i = 2:2:numstripes + for j = rowit + I = CartesianIndex(i,j) + isempty(bg.boxes[I]) || @forcefun pairwise_interactions(interaction, sim, bg.boxes[I]) + for J ∈ bg.neighborlist[I] + if !isempty(bg.boxes[J]) + @forcefun pairwise_interactions(interaction, sim, bg.boxes[I], bg.boxes[J]) + end + end # for J ∈ neighbors + end + end +end + + +@forcedef function pairwise_interactions(interaction::Function, sim::Simulation, ::NoBoxing) + @forcefun pairwise_interactions(interaction, sim, particles(sim)) +end + +@forcedef function pairwise_interactions(interaction::Function, sim::Simulation, + particles::AbstractVector{<: AbstractParticle}) + N = length(particles) + for i = 1:N + c1 = particles[i] + for j = i+1:N + @forcefun interaction(c1, particles[j], sim) + end + end +end + +@forcedef function pairwise_interactions(interaction::Function, sim::Simulation, + particles1::AbstractVector{P}, + particles2::AbstractVector{P}) where {P <: AbstractParticle} + for i in eachindex(particles1) + c1 = particles1[i] + for j in eachindex(particles2) + @forcefun interaction(c1, particles2[j], sim) + end + end +end diff --git a/src/evolve.jl b/src/evolve.jl index 71de2d25f7591b285b36a776da25bc541a44e0c0..63210135aba6667a89b44d4fc40842ef4cb65245 100644 --- a/src/evolve.jl +++ b/src/evolve.jl @@ -153,52 +153,4 @@ end computeforces!(sim) Computes the external and internal forces in the simulation. """ -@forcedef function computeforces!(sim::Simulation) - if isthreaded(sim) - @forcefun threaded_computeforces!(sim, sim.particles) - else - @forcefun serial_computeforces!(sim, sim.particles) - end -end - -# fallback for non-parallelised PCs -@forcedef threaded_computeforces!(sim, pc) = @forcefun serial_computeforces!(sim, pc) - - -""" - computeforces!(particles, sim, ::NoBoxing) -Computes all external and internal forces of `particles` naively, computing the -interactions for all pairs of particles. -""" -@forcedef function computeforces!(particles::AbstractVector{T}, - sim::Simulation, - ::NoBoxing) where {T <: AbstractParticle} - N = length(particles) - for i = 1:N - c1 = particles[i] - @forcefun internalforces!(c1, sim) - for j = i+1:N - @forcefun externalforces!(c1, particles[j], sim) - end - end -end - - - -""" - computeforces!(particles1, particles2, sim) -Computes the external forces arising due to interaction of particles from `particles1` -with particles from `particles2`. - -The function assumes that `particles1` and `particles2` do not share elements. -This condition is only checked within an [`@assert`](@ref), so there might be no -signs of things going wrong depending on the optimization level! -""" -@forcedef function computeforces!(particles1::AbstractVector{T}, particles2::AbstractVector{T}, - sim::Simulation) where {T <: AbstractParticle} - - for c1 ∈ particles1, c2 ∈ particles2 - @assert c1 != c2 - @forcefun externalforces!(c1, c2, sim) - end -end +@forcedef computeforces!(sim::Simulation) = @forcefun computeforces!(sim, sim.particles) diff --git a/src/particlecontainers/particlehandling_common.jl b/src/particlecontainers/particlehandling_common.jl index 3ad94602e0729c67f68462a528c5b399b30af405..87d33064a3bf0686ebadb1095988f7e666c1ab47 100644 --- a/src/particlecontainers/particlehandling_common.jl +++ b/src/particlecontainers/particlehandling_common.jl @@ -127,125 +127,3 @@ function findbyid(particles::AbstractVector{<:AbstractParticle}, id::UInt64) end findbyid(particles::AbstractVector{<:AbstractParticle}, ids::AbstractSet{UInt64}) = sort!([findbyid(particles, id) for id ∈ ids]) - - - -function splitbytype(particles::AbstractVector{T}) where {T<:AbstractParticle} - if !(T isa Union) - return Dict(T => particles) - end - - ret = Dict{UInt64, Vector{<:AbstractParticle}}() - for p ∈ particles - S = typeof(p) - Ssymb = hash(S) - haskey(ret, Ssymb) || (ret[Ssymb] = S[]) - push!(ret[Ssymb], p) - end - return values(ret) -end - -function maskbytype(particles::AbstractVector{U}) where {U} - if !(U isa Union) - return U, trues(length(particles)) - end - - l = length(particles) - ut = Base.uniontypes(U) - utl = length(ut) - masks = [falses(l) for i ∈ 1:utl] - - for i ∈ 1:utl - for j ∈ 1:l - masks[i][j] = isa(particles[j], ut[i]) - end - end - - return ut, masks -end - - - -## Common implementations for BoxGrid based particles containers - -""" - boxing_based_serial_computeforces(sim, pc::AbstractParticleContainer) -Uses a [`BoxGrid`](@ref) to efficiently compute pairwise interactions. -""" -@forcedef function boxing_based_serial_computeforces!(sim::Simulation, pc::AbstractParticleContainer) - # For execution order preservation compared to multithreaded simulations the - # boxes are iterated in the same striped manner. - bg = pc.boxgrid - boxparticles!(bg, particles(pc)) - - numstripes, numrows... = size(bg.boxes) - rowit = CartesianIndices(numrows) - - for i = 1:2:numstripes - for j = rowit - I = CartesianIndex(i,j) - isempty(bg.boxes[I]) || @forcefun computeforces!(bg.boxes[I], sim, NoBoxing()) - for J ∈ bg.neighborlist[I] - if !isempty(bg.boxes[J]) - @forcefun computeforces!(bg.boxes[I], bg.boxes[J], sim) - end - end # for J ∈ neighbors - end - end - for i = 2:2:numstripes - for j = rowit - I = CartesianIndex(i,j) - isempty(bg.boxes[I]) || @forcefun computeforces!(bg.boxes[I], sim, NoBoxing()) - for J ∈ bg.neighborlist[I] - if !isempty(bg.boxes[J]) - @forcefun computeforces!(bg.boxes[I], bg.boxes[J], sim) - end - end # for J ∈ neighbors - end - end -end - - -""" - boxing_based_threaded_computeforces(sim, pc::AbstractParticleContainer) -Multi-threaded version of [`boxing_based_serial_computeforces`](@ref). - -## Extended help - -Multithreading is implemented by dividing the system into two sets of interleaved stripes (slabs in 3D) -of boxes, each stripe being one box wide. -This division ensures that each stripe is independent of all other stripes in same set. - -Therefore, we can first integrate all the stripes in the first set in parallel, synchronise, then -integrate the second set. -""" -@forcedef function boxing_based_threaded_computeforces!(sim::Simulation, pc::AbstractParticleContainer) - bg = pc.boxgrid - boxparticles!(bg, particles(pc)) - - numstripes, numrows... = size(bg.boxes) - rowit = CartesianIndices(numrows) - - Threads.@threads for i = 1:2:numstripes - for j = rowit - I = CartesianIndex(i,j) - isempty(bg.boxes[I]) || @forcefun computeforces!(bg.boxes[I], sim, NoBoxing()) - for J ∈ bg.neighborlist[I] - if !isempty(bg.boxes[J]) - @forcefun computeforces!(bg.boxes[I], bg.boxes[J], sim) - end - end # for J ∈ neighbors - end - end - Threads.@threads for i = 2:2:numstripes - for j = rowit - I = CartesianIndex(i,j) - isempty(bg.boxes[I]) || @forcefun computeforces!(bg.boxes[I], sim, NoBoxing()) - for J ∈ bg.neighborlist[I] - if !isempty(bg.boxes[J]) - @forcefun computeforces!(bg.boxes[I], bg.boxes[J], sim) - end - end # for J ∈ neighbors - end - end -end diff --git a/src/particlecontainers/particleobstaclecontainer.jl b/src/particlecontainers/particleobstaclecontainer.jl index fd47eef2876606de7139bc7cc84e67478a0de431..ce31d27b7041213c4fef6b941d4a6cb157162a84 100644 --- a/src/particlecontainers/particleobstaclecontainer.jl +++ b/src/particlecontainers/particleobstaclecontainer.jl @@ -58,47 +58,30 @@ addobstacles!(ptc::ParticleObstacleContainer, obs) = append!(ptc.sov, obs) setparticles!(pc::ParticleObstacleContainer, ps) = (empty!(pc.pv); append!(pc.pv, ps)) alldynamicobjects(sim::Simulation{<:AbstractParticleContainer}) = particles(sim) - -@forcedef function serial_computeforces!(sim::Simulation, pc::ParticleObstacleContainer{N, PT, SO, BG}) where {N, PT, SO, BG<:NoBoxing} - @forcefun computeforces!(pc.pv, sim, pc.boxgrid) - for n in axes(pc.sov, 1) - so = pc.sov[n] - @forcefun obstacleforces!(so, pc.pv, sim) - end -end - - -@forcedef function serial_computeforces!(sim::Simulation, pc::ParticleObstacleContainer{N, PT, SO, BG}) where {N, PT, SO, BG<:BoxGrid} - @forcefun boxing_based_serial_computeforces!(sim, pc) - +@forcedef function computeforces!(sim::Simulation, pc::ParticleObstacleContainer) bg = pc.boxgrid - for n in axes(pc.sov, 1) - so = pc.sov[n] - for idx in pc.so_boxes[n] - @forcefun obstacleforces!(so, bg.boxes[idx], sim) - end + @forcefun pairwise_interactions(externalforces!, sim, bg) + for p in particles(sim) + internalforces!(p, sim) end -end - - -@forcedef function threaded_computeforces!(sim::Simulation, pc::ParticleObstacleContainer{N, PT, SO, BG}) where {N, PT, SO, BG<:BoxGrid} - @forcefun boxing_based_threaded_computeforces!(sim, pc) - - bg = pc.boxgrid - for n in axes(pc.sov, 1) - so = pc.sov[n] - for idx in pc.so_boxes[n] - @forcefun obstacleforces!(so, bg.boxes[idx], sim) + if pc.boxgrid isa BoxGrid + for n in axes(pc.sov, 1) + so = pc.sov[n] + for idx in pc.so_boxes[n] + # particles are sorted into boxes, within pairwise_interactions + @forcefun obstacleforces!(so, bg.boxes[idx], sim) + end + end + else + for n in axes(pc.sov, 1) + so = pc.sov[n] + @forcefun obstacleforces!(so, particles(sim), sim) end end end - - obstacles(pc::ParticleObstacleContainer) = pc.sov - - function Base.show(io::IO, ::MIME"text/plain", pc::ParticleObstacleContainer{N, PT, SO, BG}) where {N, PT, SO, BG} print(io, """ParticleObstacleContainer{$N, $PT, $SO, $BG} diff --git a/src/particlecontainers/polocontainer.jl b/src/particlecontainers/polocontainer.jl index 44202a2c775f3283c568d5cd2f361aef4cd5aeef..aa6d7c4fcf85780c0edf15dc076bbd95e80fee30 100644 --- a/src/particlecontainers/polocontainer.jl +++ b/src/particlecontainers/polocontainer.jl @@ -60,33 +60,17 @@ setparticles!(pc::POLOContainer, ps) = (empty!(pc.pv); append!(pc.pv, ps)) alldynamicobjects(sim::Simulation{<:POLOContainer}) = vcat(sim.particles.pv, sim.particles.lov) -@forcedef function serial_computeforces!(sim::Simulation, pc::POLOContainer) - @forcefun boxing_based_serial_computeforces!(sim, pc) - - bg = pc.boxgrid - for n in axes(pc.sov, 1) - so = pc.sov[n] - for idx in pc.so_boxes[n] - @forcefun obstacleforces!(so, bg.boxes[idx], sim) - end - end - - for n in axes(pc.lov, 1) - lo = pc.lov[n] - for cell in particles(sim) - @forcefun externalforces!(lo, cell, sim) - end +@forcedef function computeforces!(sim::Simulation, pc::POLOContainer) + @forcefun pairwise_interactions(externalforces!, sim, pc.boxgrid) + for p in particles(sim) + internalforces!(p, sim) end -end - - -@forcedef function threaded_computeforces!(sim::Simulation, pc::POLOContainer) - @forcefun boxing_based_threaded_computeforces!(sim, pc) bg = pc.boxgrid for n in axes(pc.sov, 1) so = pc.sov[n] for idx in pc.so_boxes[n] + # particles are sorted into boxes, within pairwise_interactions @forcefun obstacleforces!(so, bg.boxes[idx], sim) end end @@ -99,7 +83,6 @@ end end end - obstacles(pc::POLOContainer) = pc.sov diff --git a/src/particlecontainers/simpleparticlecontainer.jl b/src/particlecontainers/simpleparticlecontainer.jl index d582f6e200b9e46564aea215fba1d52a15761c92..cd66c9870b21ce5e329ff00d6b5948f4a27182c6 100644 --- a/src/particlecontainers/simpleparticlecontainer.jl +++ b/src/particlecontainers/simpleparticlecontainer.jl @@ -47,22 +47,13 @@ alldynamicobjects(sim::Simulation{<:SimpleParticleContainer}) = sim.particles.pa ## FORCE COMPUTATION ########################################################### -@forcedef function serial_computeforces!(sim::Simulation, pc::SimpleParticleContainer{N, PT, BG}) where {N, PT, BG<:NoBoxing} - @forcefun computeforces!(pc.particles, sim, pc.boxgrid) +@forcedef function computeforces!(sim::Simulation, pc::SimpleParticleContainer) + @forcefun pairwise_interactions(externalforces!, sim, pc.boxgrid) + for p in particles(sim) + internalforces!(p, sim) + end end - -@forcedef function serial_computeforces!(sim::Simulation, pc::SimpleParticleContainer{N, PT, BG}) where {N, PT, BG<:BoxGrid} - @forcefun boxing_based_serial_computeforces!(sim, pc) -end - -@forcedef function threaded_computeforces!(sim::Simulation, pc::SimpleParticleContainer{N, PT, BG}) where {N, PT, BG<:BoxGrid} - @forcefun boxing_based_threaded_computeforces!(sim, pc) -end - - - - function Base.show(io::IO, ::MIME"text/plain", pc::SimpleParticleContainer{N, PT, BG}) where {N, PT, BG} print(io, """SimpleParticleContainer{$N, $PT, $BG}