From 129f38b6433316416acd8f7d88c0ec5b68e2b38e Mon Sep 17 00:00:00 2001 From: Jonas Isensee <jonas.isensee@ds.mpg.de> Date: Mon, 18 Sep 2023 11:03:13 +0200 Subject: [PATCH 1/2] new pairwise_interaction function --- src/boxgrid.jl | 105 +++++++++++++++ src/evolve.jl | 50 +------ .../particlehandling_common.jl | 122 ------------------ .../particleobstaclecontainer.jl | 49 +++---- src/particlecontainers/polocontainer.jl | 27 +--- .../simpleparticlecontainer.jl | 19 +-- 6 files changed, 132 insertions(+), 240 deletions(-) diff --git a/src/boxgrid.jl b/src/boxgrid.jl index 8118365..ed2055e 100644 --- a/src/boxgrid.jl +++ b/src/boxgrid.jl @@ -158,3 +158,108 @@ 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 + +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. +""" +function pairwise_interactions(interaction::Functions, 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, bg.boxes[I], sim) + for J ∈ bg.neighborlist[I] + if !isempty(bg.boxes[J]) + @forcefun pairwise_interactions(interaction, 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 pairwise_interactions(interaction, bg.boxes[I], sim) + for J ∈ bg.neighborlist[I] + if !isempty(bg.boxes[J]) + @forcefun pairwise_interactions(interaction, bg.boxes[I], bg.boxes[J], sim) + 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::Functions, 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 diff --git a/src/evolve.jl b/src/evolve.jl index 71de2d2..6321013 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 3ad9460..87d3306 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 fd47eef..ce31d27 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 44202a2..aa6d7c4 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 d582f6e..cd66c98 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} -- GitLab From 50100fbbe73afc1de989bd17dba1cf38859f3b85 Mon Sep 17 00:00:00 2001 From: Jonas Isensee <jonas.isensee@ds.mpg.de> Date: Mon, 18 Sep 2023 13:14:48 +0200 Subject: [PATCH 2/2] fix definitions --- src/boxgrid.jl | 25 ++++++++++++++++++------- 1 file changed, 18 insertions(+), 7 deletions(-) diff --git a/src/boxgrid.jl b/src/boxgrid.jl index ed2055e..ec9cfa9 100644 --- a/src/boxgrid.jl +++ b/src/boxgrid.jl @@ -163,7 +163,7 @@ end ## BoxGrid pairwise interaction computation -pairwise_interactions(interaction::Function, sim::Simulation, bg::BoxGrid) = +@forcedef pairwise_interactions(interaction::Function, sim::Simulation, bg::BoxGrid) = pairwise_interactions(interaction, sim, bg, Val{isthreaded(sim)}()) @@ -172,7 +172,7 @@ pairwise_interactions(interaction::Function, sim::Simulation, bg::BoxGrid) = Uses a [`BoxGrid`](@ref) to efficiently compute pairwise interactions. `BoxGrid` structure is mutated in the sense that particles are sorted into their boxes. """ -function pairwise_interactions(interaction::Functions, sim::Simulation, bg::BoxGrid, ::Val{false}) +@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)) @@ -183,10 +183,10 @@ function pairwise_interactions(interaction::Functions, sim::Simulation, bg::BoxG for i = 1:2:numstripes for j = rowit I = CartesianIndex(i,j) - isempty(bg.boxes[I]) || @forcefun pairwise_interactions(interaction, bg.boxes[I], sim) + 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, bg.boxes[I], bg.boxes[J], sim) + @forcefun pairwise_interactions(interaction, sim, bg.boxes[I], bg.boxes[J]) end end # for J ∈ neighbors end @@ -194,10 +194,10 @@ function pairwise_interactions(interaction::Functions, sim::Simulation, bg::BoxG for i = 2:2:numstripes for j = rowit I = CartesianIndex(i,j) - isempty(bg.boxes[I]) || @forcefun pairwise_interactions(interaction, bg.boxes[I], sim) + 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, bg.boxes[I], bg.boxes[J], sim) + @forcefun pairwise_interactions(interaction, sim, bg.boxes[I], bg.boxes[J]) end end # for J ∈ neighbors end @@ -218,7 +218,7 @@ This division ensures that each stripe is independent of all other stripes in sa 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::Functions, sim::Simulation, bg::BoxGrid, ::Val{true}) +@forcedef function pairwise_interactions(interaction::Function, sim::Simulation, bg::BoxGrid, ::Val{true}) boxparticles!(bg, particles(sim)) numstripes, numrows... = size(bg.boxes) @@ -263,3 +263,14 @@ end 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 -- GitLab