diff --git a/src/SPHCellList.jl b/src/SPHCellList.jl index 19ddae69..760678e4 100644 --- a/src/SPHCellList.jl +++ b/src/SPHCellList.jl @@ -18,7 +18,10 @@ using ..OpenExternalPrograms using ..SPHKernels using ..SPHViscosityModels using ..SPHDensityDiffusionModels -using ..SPHNeighborList: BuildNeighborCellLists!, ComputeCellNeighborCounts, ComputeCellParticleCounts, ConstructStencil, ExtractCells!, MapFloor, UpdateNeighbors!, UpdateΔx! +using ..SPHNeighborList: BuildNeighborCellLists!, CellLookup, DenseCellLookup, ComputeCellNeighborCounts!, + ComputeCellParticleCounts!, ConstructStencil, ExtractCells!, GetCellIndexDispatch, + InitializeCellLookup, InitializeDenseCellLookup, MapFloor, NeighborListScratch, + UpdateCellCounts!, UpdateNeighbors!, UpdateΔx! using StaticArrays using StructArrays: StructArray, foreachfield @@ -35,7 +38,7 @@ using LinearAlgebra function NeighborLoopPerParticle!(SimDensityDiffusion::SDD, SimViscosity::SV, SimKernel, SimMetaData::SimulationMetaData{D,T,NoShifting,NoKernelOutput,B,L}, SimConstants, SimParticles, ParticleRanges, - CellDict, NeighborCellLists, dρdtI, + CellIndices, NeighborCellOffsets, NeighborCellCounts, NeighborCells, dρdtI, Acceleration, ∇Cᵢ, ∇◌rᵢ, max_visc = nothing, min_dt_force = nothing; @@ -51,12 +54,9 @@ using LinearAlgebra @inbounds Threads.@threads for i in eachindex(Position) dρdt_acc = zero(dρdtI[i]) acc_acc = zero(Acceleration[i]) - CellIndex = Cells[i] - CellListIndex = get(CellDict, CellIndex, 1) + CellListIndex = CellIndices[i] SameCellStart = ParticleRanges[CellListIndex] SameCellEnd = ParticleRanges[CellListIndex + 1] - 1 - NeighborCellIndices = NeighborCellLists[CellListIndex] - @inbounds for j in SameCellStart:SameCellEnd jIndex = ParticleOrder[j] if jIndex != i @@ -67,7 +67,10 @@ using LinearAlgebra ) end end - for NeighborIdx in NeighborCellIndices + neighbor_start = NeighborCellOffsets[CellListIndex] + neighbor_count = NeighborCellCounts[CellListIndex] + for n in 0:(neighbor_count - 1) + NeighborIdx = NeighborCells[neighbor_start + n] StartIndex_ = ParticleRanges[NeighborIdx] EndIndex_ = ParticleRanges[NeighborIdx + 1] - 1 @inbounds for j in StartIndex_:EndIndex_ @@ -91,7 +94,7 @@ using LinearAlgebra function NeighborLoopPerParticle!(SimDensityDiffusion::SDD, SimViscosity::SV, SimKernel, SimMetaData::SimulationMetaData{D,T,NoShifting,K,B,L}, SimConstants, SimParticles, ParticleRanges, - CellDict, NeighborCellLists, dρdtI, + CellIndices, NeighborCellOffsets, NeighborCellCounts, NeighborCells, dρdtI, Acceleration, ∇Cᵢ, ∇◌rᵢ, max_visc = nothing, min_dt_force = nothing; @@ -110,12 +113,9 @@ using LinearAlgebra acc_acc = zero(Acceleration[i]) kernel_acc = zero(Kernel[i]) kernel_grad_acc = zero(KernelGradient[i]) - CellIndex = Cells[i] - CellListIndex = get(CellDict, CellIndex, 1) + CellListIndex = CellIndices[i] SameCellStart = ParticleRanges[CellListIndex] SameCellEnd = ParticleRanges[CellListIndex + 1] - 1 - NeighborCellIndices = NeighborCellLists[CellListIndex] - @inbounds for j in SameCellStart:SameCellEnd jIndex = ParticleOrder[j] if jIndex != i @@ -128,7 +128,10 @@ using LinearAlgebra ) end end - for NeighborIdx in NeighborCellIndices + neighbor_start = NeighborCellOffsets[CellListIndex] + neighbor_count = NeighborCellCounts[CellListIndex] + for n in 0:(neighbor_count - 1) + NeighborIdx = NeighborCells[neighbor_start + n] StartIndex_ = ParticleRanges[NeighborIdx] EndIndex_ = ParticleRanges[NeighborIdx + 1] - 1 @inbounds for j in StartIndex_:EndIndex_ @@ -156,7 +159,7 @@ using LinearAlgebra function NeighborLoopPerParticle!(SimDensityDiffusion::SDD, SimViscosity::SV, SimKernel, SimMetaData::SimulationMetaData{D,T,S,NoKernelOutput,B,L}, SimConstants, SimParticles, ParticleRanges, - CellDict, NeighborCellLists, dρdtI, + CellIndices, NeighborCellOffsets, NeighborCellCounts, NeighborCells, dρdtI, Acceleration, ∇Cᵢ, ∇◌rᵢ, max_visc = nothing, min_dt_force = nothing; @@ -174,12 +177,9 @@ using LinearAlgebra acc_acc = zero(Acceleration[i]) shift_c_acc = zero(∇Cᵢ[i]) shift_r_acc = zero(∇◌rᵢ[i]) - CellIndex = Cells[i] - CellListIndex = get(CellDict, CellIndex, 1) + CellListIndex = CellIndices[i] SameCellStart = ParticleRanges[CellListIndex] SameCellEnd = ParticleRanges[CellListIndex + 1] - 1 - NeighborCellIndices = NeighborCellLists[CellListIndex] - @inbounds for j in SameCellStart:SameCellEnd jIndex = ParticleOrder[j] if jIndex != i @@ -192,7 +192,10 @@ using LinearAlgebra ) end end - for NeighborIdx in NeighborCellIndices + neighbor_start = NeighborCellOffsets[CellListIndex] + neighbor_count = NeighborCellCounts[CellListIndex] + for n in 0:(neighbor_count - 1) + NeighborIdx = NeighborCells[neighbor_start + n] StartIndex_ = ParticleRanges[NeighborIdx] EndIndex_ = ParticleRanges[NeighborIdx + 1] - 1 @inbounds for j in StartIndex_:EndIndex_ @@ -220,7 +223,7 @@ using LinearAlgebra function NeighborLoopPerParticle!(SimDensityDiffusion::SDD, SimViscosity::SV, SimKernel, SimMetaData::SimulationMetaData{D,T,S,K,B,L}, SimConstants, SimParticles, ParticleRanges, - CellDict, NeighborCellLists, dρdtI, + CellIndices, NeighborCellOffsets, NeighborCellCounts, NeighborCells, dρdtI, Acceleration, ∇Cᵢ, ∇◌rᵢ, max_visc = nothing, min_dt_force = nothing; @@ -242,12 +245,9 @@ using LinearAlgebra kernel_grad_acc = zero(KernelGradient[i]) shift_c_acc = zero(∇Cᵢ[i]) shift_r_acc = zero(∇◌rᵢ[i]) - CellIndex = Cells[i] - CellListIndex = get(CellDict, CellIndex, 1) + CellListIndex = CellIndices[i] SameCellStart = ParticleRanges[CellListIndex] SameCellEnd = ParticleRanges[CellListIndex + 1] - 1 - NeighborCellIndices = NeighborCellLists[CellListIndex] - @inbounds for j in SameCellStart:SameCellEnd jIndex = ParticleOrder[j] if jIndex != i @@ -260,7 +260,10 @@ using LinearAlgebra ) end end - for NeighborIdx in NeighborCellIndices + neighbor_start = NeighborCellOffsets[CellListIndex] + neighbor_count = NeighborCellCounts[CellListIndex] + for n in 0:(neighbor_count - 1) + NeighborIdx = NeighborCells[neighbor_start + n] StartIndex_ = ParticleRanges[NeighborIdx] EndIndex_ = ParticleRanges[NeighborIdx + 1] - 1 @inbounds for j in StartIndex_:EndIndex_ @@ -291,7 +294,7 @@ using LinearAlgebra f(SimKernel, GhostPoint) = CartesianIndex(map(x -> MapFloor(x, SimKernel.H⁻¹), Tuple(GhostPoint))) function NeighborLoopMDBC!(SimKernel, SimMetaData::SimulationMetaData{Dimensions, FloatType, SMode, KMode, BMode, LMode}, - SimConstants, ParticleRanges, CellDict, Position, + SimConstants, ParticleRanges, CellLookup, DenseLookup, UseDense, Position, Density, GhostPoints, GhostNormals, ParticleType, bᵧ, Aᵧ; ParticleOrder) where {Dimensions, FloatType, SMode, KMode, BMode, LMode} @@ -314,7 +317,7 @@ using LinearAlgebra # Returns a range, x>:x for exact match and x=:x for no match # utilizes that it is a sorted array and requires no isequal constructor, # so I prefer this for now - NeighborIdx = get(CellDict, SCellIndex, 1) + NeighborIdx = GetCellIndexDispatch(UseDense, DenseLookup, CellLookup, SCellIndex, 1) StartIndex_ = ParticleRanges[NeighborIdx] EndIndex_ = ParticleRanges[NeighborIdx + 1] - 1 @@ -632,7 +635,7 @@ using LinearAlgebra function ApplyMDBCBeforeHalf!(SimMetaData::SimulationMetaData{D,T,S,K,SimpleMDBC,L}, SimKernel, SimConstants, SimParticles, - ParticleRanges, CellDict, Position, Density, + ParticleRanges, CellLookup, DenseLookup, UseDense, Position, Density, GhostPoints, GhostNormals, ParticleType; ParticleOrder ) where {D,T,S<:ShiftingMode,K<:KernelOutputMode,L<:LogMode} @@ -640,7 +643,7 @@ using LinearAlgebra DimensionsPlus = D + 1 bᵧ = @alloc(SVector{DimensionsPlus, T}, length(Position)) Aᵧ = @alloc(SMatrix{DimensionsPlus, DimensionsPlus, T, DimensionsPlus*DimensionsPlus}, length(Position)) - NeighborLoopMDBC!(SimKernel, SimMetaData, SimConstants, ParticleRanges, CellDict, Position, Density, GhostPoints, GhostNormals, ParticleType, bᵧ, Aᵧ; ParticleOrder = ParticleOrder) + NeighborLoopMDBC!(SimKernel, SimMetaData, SimConstants, ParticleRanges, CellLookup, DenseLookup, UseDense, Position, Density, GhostPoints, GhostNormals, ParticleType, bᵧ, Aᵧ; ParticleOrder = ParticleOrder) ApplyMDBCCorrection(SimConstants, SimParticles, bᵧ, Aᵧ) end @@ -697,10 +700,10 @@ using LinearAlgebra @inbounds function SimulationLoop(SimDensityDiffusion::SDD, SimViscosity::SV, SimKernel, SimMetaData::SimulationMetaData{Dimensions, FloatType, SMode, KMode, BMode, LMode}, SimConstants, SimParticles, FullStencil, - ParticleRanges, UniqueCells, CellDict, - ParticleOrder, CellOffsets, - NeighborCellLists, dρdtI, Velocityₙ⁺, - Positionₙ⁺, ρₙ⁺, ∇Cᵢ, ∇◌rᵢ, + ParticleRanges, UniqueCells, CellLookup, DenseLookup, UseDense, + ParticleOrder, CellOffsets, CellIndices, + NeighborCellOffsets, NeighborCellCounts, NeighborCells, + NeighborScratch, dρdtI, Velocityₙ⁺, Positionₙ⁺, ρₙ⁺, ∇Cᵢ, ∇◌rᵢ, MotionDefinition::Union{ Nothing, AbstractVector{ @@ -747,10 +750,34 @@ using LinearAlgebra # Remove if statement logic if you want to update each iteration # if mod(SimMetaData.Iteration, ceil(Int, SimKernel.H / (SimConstants.c₀ * dt * (1/SimConstants.CFL)) )) == 0 || SimMetaData.Iteration == 1 if ShouldRebuild - @timeit SimMetaData.HourGlass "01a Actual Calculate IndexCounter" SimMetaData.IndexCounter = UpdateNeighbors!(SimParticles, SimKernel.H⁻¹, ParticleRanges, UniqueCells, CellDict, ParticleOrder, CellOffsets) + @timeit SimMetaData.HourGlass "01a Actual Calculate IndexCounter" begin + SimMetaData.IndexCounter, UseDense[] = UpdateNeighbors!( + SimParticles, + SimKernel.H⁻¹, + ParticleRanges, + UniqueCells, + CellLookup, + ParticleOrder, + CellOffsets, + CellIndices, + DenseLookup, + ) + end + UpdateCellCounts!(NeighborScratch, ParticleRanges, SimMetaData.IndexCounter) SimMetaData.Δx = zero(eltype(dρdtI)) UniqueCellsView = view(UniqueCells, 1:SimMetaData.IndexCounter) - BuildNeighborCellLists!(NeighborCellLists, FullStencil, UniqueCellsView, ParticleRanges, CellDict) + BuildNeighborCellLists!( + NeighborCellOffsets, + NeighborCellCounts, + NeighborCells, + FullStencil, + UniqueCellsView, + ParticleRanges, + CellLookup, + DenseLookup, + Val(UseDense[]), + NeighborScratch.CellOccupied, + ) end end @@ -760,21 +787,23 @@ using LinearAlgebra if SimMetaData isa SimulationMetaData{Dimensions, FloatType, SMode, KMode, NoMDBC, LMode} where {SMode, KMode, LMode} @timeit SimMetaData.HourGlass "03 Apply MDBC before Half TimeStep" ApplyMDBCBeforeHalf!(SimMetaData) else - @timeit SimMetaData.HourGlass "03 Apply MDBC before Half TimeStep" ApplyMDBCBeforeHalf!(SimMetaData, SimKernel, SimConstants, SimParticles, ParticleRanges, CellDict, Position, Density, GhostPoints, GhostNormals, ParticleType; ParticleOrder = ParticleOrder) + @timeit SimMetaData.HourGlass "03 Apply MDBC before Half TimeStep" ApplyMDBCBeforeHalf!(SimMetaData, SimKernel, SimConstants, SimParticles, ParticleRanges, CellLookup, DenseLookup, Val(UseDense[]), Position, Density, GhostPoints, GhostNormals, ParticleType; ParticleOrder = ParticleOrder) end if SimMetaData isa SimulationMetaData{Dimensions, FloatType, SMode, NoKernelOutput, BMode, LMode} where {SMode, BMode, LMode} @timeit SimMetaData.HourGlass "04 First NeighborLoop" NeighborLoopPerParticle!( SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData, - SimConstants, SimParticles, ParticleRanges, CellDict, - NeighborCellLists, dρdtI, Acceleration, ∇Cᵢ, ∇◌rᵢ, + SimConstants, SimParticles, ParticleRanges, CellIndices, + NeighborCellOffsets, NeighborCellCounts, NeighborCells, dρdtI, + Acceleration, ∇Cᵢ, ∇◌rᵢ, ParticleOrder = ParticleOrder, ) else @timeit SimMetaData.HourGlass "04 First NeighborLoop" NeighborLoopPerParticle!( SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData, - SimConstants, SimParticles, ParticleRanges, CellDict, - NeighborCellLists, dρdtI, Acceleration, ∇Cᵢ, ∇◌rᵢ, + SimConstants, SimParticles, ParticleRanges, CellIndices, + NeighborCellOffsets, NeighborCellCounts, NeighborCells, dρdtI, + Acceleration, ∇Cᵢ, ∇◌rᵢ, ParticleOrder = ParticleOrder, ) end @@ -791,8 +820,9 @@ using LinearAlgebra if SimMetaData isa SimulationMetaData{Dimensions, FloatType, SMode, NoKernelOutput, BMode, LMode} where {SMode, BMode, LMode} @timeit SimMetaData.HourGlass "08 Second NeighborLoop" NeighborLoopPerParticle!( SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData, - SimConstants, SimParticles, ParticleRanges, CellDict, - NeighborCellLists, dρdtI, Acceleration, ∇Cᵢ, ∇◌rᵢ, + SimConstants, SimParticles, ParticleRanges, CellIndices, + NeighborCellOffsets, NeighborCellCounts, NeighborCells, dρdtI, + Acceleration, ∇Cᵢ, ∇◌rᵢ, max_visc, min_dt_force, ParticleOrder = ParticleOrder, Position = Positionₙ⁺, @@ -802,8 +832,9 @@ using LinearAlgebra else @timeit SimMetaData.HourGlass "08 Second NeighborLoop" NeighborLoopPerParticle!( SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData, - SimConstants, SimParticles, ParticleRanges, CellDict, - NeighborCellLists, dρdtI, Acceleration, ∇Cᵢ, ∇◌rᵢ, + SimConstants, SimParticles, ParticleRanges, CellIndices, + NeighborCellOffsets, NeighborCellCounts, NeighborCells, dρdtI, + Acceleration, ∇Cᵢ, ∇◌rᵢ, max_visc, min_dt_force, ParticleOrder = ParticleOrder, Position = Positionₙ⁺, @@ -828,17 +859,35 @@ using LinearAlgebra end # To control if grid is exported in output - function DetermineOutput(::Val{true}, SimMetaData, output, ParticleRanges, NeighborCellLists, UniqueCellsView) - cell_particle_counts = ComputeCellParticleCounts(ParticleRanges, length(UniqueCellsView)) - cell_neighbor_counts = ComputeCellNeighborCounts(ParticleRanges, NeighborCellLists, length(UniqueCellsView)) + function DetermineOutput(::Val{true}, SimMetaData, output, ParticleRanges, + NeighborCellOffsets, NeighborCellCounts, NeighborCells, + NeighborScratch, + UniqueCellsView) + UpdateCellCounts!(NeighborScratch, ParticleRanges, length(UniqueCellsView)) + ComputeCellNeighborCounts!( + NeighborScratch.NeighborCounts, + NeighborScratch.CellCounts, + NeighborCellOffsets, + NeighborCellCounts, + NeighborCells, + length(UniqueCellsView), + ) @timeit SimMetaData.HourGlass "13 Save Particle Data" begin output.enqueue_particles(SimMetaData.OutputIterationCounter) - output.enqueue_grid(SimMetaData.OutputIterationCounter, UniqueCellsView, cell_particle_counts=cell_particle_counts, cell_neighbor_counts=cell_neighbor_counts) + output.enqueue_grid( + SimMetaData.OutputIterationCounter, + UniqueCellsView, + cell_particle_counts=NeighborScratch.CellCounts, + cell_neighbor_counts=NeighborScratch.NeighborCounts, + ) end end - function DetermineOutput(::Val{false}, SimMetaData, output, ParticleRanges, NeighborCellLists, UniqueCellsView) + function DetermineOutput(::Val{false}, SimMetaData, output, ParticleRanges, + NeighborCellOffsets, NeighborCellCounts, NeighborCells, + NeighborScratch, + UniqueCellsView) output.enqueue_particles(SimMetaData.OutputIterationCounter) end @@ -867,11 +916,17 @@ using LinearAlgebra # Produce sorting related variables ParticleRanges = zeros(Int, NumberOfPoints + 1 + 1) # +1 for the last particle, +1 for dummy entry UniqueCells = zeros(CartesianIndex{Dimensions}, NumberOfPoints) - CellDict = Dict{CartesianIndex{Dimensions}, Int}() + CellLookup = InitializeCellLookup(Val(Dimensions), NumberOfPoints * 2) + DenseLookup = InitializeDenseCellLookup(Val(Dimensions)) + UseDense = Ref(false) FullStencil = ConstructStencil(Val(Dimensions)) - NeighborCellLists = [Int[] for _ in 1:length(UniqueCells)] + NeighborCellOffsets = zeros(Int, length(UniqueCells)) + NeighborCellCounts = zeros(Int, length(UniqueCells)) + NeighborCells = zeros(Int, length(UniqueCells) * (length(FullStencil) - 1)) + NeighborScratch = NeighborListScratch() ParticleOrder = zeros(Int, NumberOfPoints) CellOffsets = zeros(Int, length(ParticleRanges)) + CellIndices = zeros(Int, NumberOfPoints) output = SetupVTKOutput(SimMetaData, SimParticles, SimKernel, Dimensions) @@ -883,15 +938,17 @@ using LinearAlgebra cell_particle_counts = nothing cell_neighbor_counts = nothing if SimMetaData.ExportGridCellParticleCounts - cell_particle_counts = ComputeCellParticleCounts( - ParticleRanges, - SimMetaData.IndexCounter, - ) - cell_neighbor_counts = ComputeCellNeighborCounts( - ParticleRanges, - NeighborCellLists, + UpdateCellCounts!(NeighborScratch, ParticleRanges, SimMetaData.IndexCounter) + ComputeCellNeighborCounts!( + NeighborScratch.NeighborCounts, + NeighborScratch.CellCounts, + NeighborCellOffsets, + NeighborCellCounts, + NeighborCells, SimMetaData.IndexCounter, ) + cell_particle_counts = NeighborScratch.CellCounts + cell_neighbor_counts = NeighborScratch.NeighborCounts end output.enqueue_grid( SimMetaData.OutputIterationCounter, @@ -909,8 +966,9 @@ using LinearAlgebra @timeit SimMetaData.HourGlass "00 SimulationLoop" SimulationLoop( SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData, SimConstants, SimParticles, FullStencil, ParticleRanges, - UniqueCells, CellDict, ParticleOrder, CellOffsets, - NeighborCellLists, dρdtI, Velocityₙ⁺, Positionₙ⁺, ρₙ⁺, + UniqueCells, CellLookup, DenseLookup, UseDense, ParticleOrder, CellOffsets, CellIndices, + NeighborCellOffsets, NeighborCellCounts, NeighborCells, NeighborScratch, + dρdtI, Velocityₙ⁺, Positionₙ⁺, ρₙ⁺, ∇Cᵢ, ∇◌rᵢ, MotionDefinition, ) push!(SimMetaData.TimeSteps, SimMetaData.CurrentTimeStep) @@ -921,7 +979,17 @@ using LinearAlgebra UniqueCellsView = view(UniqueCells, 1:SimMetaData.IndexCounter) - @timeit SimMetaData.HourGlass "13 Determine Output" DetermineOutput(Val(SimMetaData.ExportGridCellParticleCounts), SimMetaData,output, ParticleRanges, NeighborCellLists, UniqueCellsView) + @timeit SimMetaData.HourGlass "13 Determine Output" DetermineOutput( + Val(SimMetaData.ExportGridCellParticleCounts), + SimMetaData, + output, + ParticleRanges, + NeighborCellOffsets, + NeighborCellCounts, + NeighborCells, + NeighborScratch, + UniqueCellsView, + ) if SimMetaData.TotalTime > SimMetaData.SimulationTime diff --git a/src/SPHExample.jl b/src/SPHExample.jl index a56895ea..98bb62e6 100644 --- a/src/SPHExample.jl +++ b/src/SPHExample.jl @@ -58,7 +58,9 @@ module SPHExample using .SPHNeighborList export ConstructStencil, ExtractCells!, UpdateNeighbors!, - BuildNeighborCellLists!, ComputeCellParticleCounts, ComputeCellNeighborCounts + BuildNeighborCellLists!, NeighborListScratch, UpdateCellCounts!, + ComputeCellParticleCounts!, ComputeCellNeighborCounts!, + ComputeCellParticleCounts, ComputeCellNeighborCounts using .SPHCellList export NeighborLoop!, ComputeInteractions!, RunSimulation diff --git a/src/SPHNeighborList.jl b/src/SPHNeighborList.jl index b94fc07e..64bce3ff 100644 --- a/src/SPHNeighborList.jl +++ b/src/SPHNeighborList.jl @@ -1,41 +1,212 @@ module SPHNeighborList export ConstructStencil, ExtractCells!, UpdateNeighbors!, - BuildNeighborCellLists!, ComputeCellParticleCounts, ComputeCellNeighborCounts, UpdateΔx! + BuildNeighborCellLists!, ComputeCellParticleCounts, ComputeCellNeighborCounts, UpdateΔx!, + CellLookup, DenseCellLookup, InitializeCellLookup, InitializeDenseCellLookup, + ResetCellLookup!, GetCellIndex, SetCellIndex!, GetCellIndexDense, GetCellIndexDispatch, + NeighborListScratch, UpdateCellCounts!, ComputeCellParticleCounts!, ComputeCellNeighborCounts! using StaticArrays +mutable struct CellLookup{D} + Keys::Vector{CartesianIndex{D}} + Values::Vector{Int} + Stamps::Vector{UInt32} + Mask::Int + Epoch::UInt32 +end + +mutable struct DenseCellLookup{D} + MinCell::SVector{D, Int} + Dims::SVector{D, Int} + Strides::SVector{D, Int} + IndexMap::Vector{Int} +end + +mutable struct NeighborListScratch + CellCounts::Vector{Int} + CellOccupied::Vector{Bool} + NeighborCounts::Vector{Int} + + function NeighborListScratch() + new(Int[], Bool[], Int[]) + end +end + +@inline function ResizeNeighborListScratch!(Scratch::NeighborListScratch, CellCount) + resize!(Scratch.CellCounts, CellCount) + resize!(Scratch.CellOccupied, CellCount) + resize!(Scratch.NeighborCounts, CellCount) + return nothing +end + +@inline function ZeroCell(::Val{D}) where D + return CartesianIndex(ntuple(_ -> 0, Val(D))) +end + +@inline function NextPow2(n::Int) + size = 1 + while size < n + size <<= 1 + end + return size +end + +function InitializeCellLookup(::Val{D}, capacity::Int) where D + size = NextPow2(max(16, capacity)) + zero_cell = ZeroCell(Val(D)) + return CellLookup{D}(fill(zero_cell, size), zeros(Int, size), zeros(UInt32, size), size - 1, UInt32(1)) +end + +function InitializeDenseCellLookup(::Val{D}) where D + zero_vec = SVector{D, Int}(ntuple(_ -> 0, Val(D))) + return DenseCellLookup{D}(zero_vec, zero_vec, zero_vec, Int[]) +end + +@inline function DenseIndex(Dims, Strides, MinCell, Cell) + index = 1 + @inbounds for i in 1:length(Dims) + index += (Cell[i] - MinCell[i]) * Strides[i] + end + return index +end + +@inline function GetCellIndexDense(Lookup::DenseCellLookup{D}, Cell::CartesianIndex{D}, default::Int) where D + MinCell = Lookup.MinCell + Dims = Lookup.Dims + @inbounds for i in 1:D + v = Cell.I[i] - MinCell[i] + if v < 0 || v >= Dims[i] + return default + end + end + index = DenseIndex(Dims, Lookup.Strides, MinCell, Cell.I) + value = Lookup.IndexMap[index] + return value == 0 ? default : value +end + +@inline function GetCellIndexDispatch(::Val{true}, DenseLookup, CellLookup, Cell, default) + return GetCellIndexDense(DenseLookup, Cell, default) +end + +@inline function GetCellIndexDispatch(::Val{false}, DenseLookup, CellLookup, Cell, default) + return GetCellIndex(CellLookup, Cell, default) +end + +function ResetCellLookup!(Lookup::CellLookup{D}, capacity::Int) where D + target = max(16, capacity) + if length(Lookup.Keys) < target + size = NextPow2(target) + zero_cell = ZeroCell(Val(D)) + Lookup.Keys = fill(zero_cell, size) + Lookup.Values = zeros(Int, size) + Lookup.Stamps = zeros(UInt32, size) + Lookup.Mask = size - 1 + Lookup.Epoch = UInt32(1) + else + Lookup.Epoch += 1 + if Lookup.Epoch == UInt32(0) + fill!(Lookup.Stamps, 0) + Lookup.Epoch = UInt32(1) + end + end + return nothing +end + +@inline function HashCellIndex(Cell::CartesianIndex{D}) where D + h = UInt(0x9e3779b97f4a7c15) + @inbounds for i in 1:D + v = unsigned(Cell.I[i]) + h ⊻= v + 0x9e3779b97f4a7c15 + (h << 6) + (h >> 2) + end + return h +end + +@inline function GetCellIndex(Lookup::CellLookup{D}, Cell::CartesianIndex{D}, default::Int) where D + mask = Lookup.Mask + index = Int(HashCellIndex(Cell) & UInt(mask)) + @inbounds for _ in 0:mask + slot = index + 1 + if Lookup.Stamps[slot] != Lookup.Epoch + return default + elseif Lookup.Keys[slot] == Cell + return Lookup.Values[slot] + end + index = (index + 1) & mask + end + return default +end + +@inline function SetCellIndex!(Lookup::CellLookup{D}, Cell::CartesianIndex{D}, value::Int) where D + mask = Lookup.Mask + index = Int(HashCellIndex(Cell) & UInt(mask)) + @inbounds for _ in 0:mask + slot = index + 1 + if Lookup.Stamps[slot] != Lookup.Epoch + Lookup.Stamps[slot] = Lookup.Epoch + Lookup.Keys[slot] = Cell + Lookup.Values[slot] = value + return nothing + elseif Lookup.Keys[slot] == Cell + Lookup.Values[slot] = value + return nothing + end + index = (index + 1) & mask + end + return nothing +end + +@inline function GetOrInsertCellIndex!(Lookup::CellLookup{D}, + Cell::CartesianIndex{D}, + NextIndex::Int) where D + mask = Lookup.Mask + index = Int(HashCellIndex(Cell) & UInt(mask)) + @inbounds for _ in 0:mask + slot = index + 1 + if Lookup.Stamps[slot] != Lookup.Epoch + Lookup.Stamps[slot] = Lookup.Epoch + Lookup.Keys[slot] = Cell + Lookup.Values[slot] = NextIndex + return NextIndex, true + elseif Lookup.Keys[slot] == Cell + return Lookup.Values[slot], false + end + index = (index + 1) & mask + end + return NextIndex, true +end + function ConstructStencil(V::Val{d}) where d return CartesianIndices(ntuple(_ -> -1:1, V)) end -function BuildNeighborCellLists!(NeighborCellLists, FullStencil, UniqueCellsView, ParticleRanges, CellDict) - TargetLen = length(UniqueCellsView) - OriginalLen = length(NeighborCellLists) - resize!(NeighborCellLists, TargetLen) +function BuildNeighborCellLists!(NeighborCellOffsets, NeighborCellCounts, NeighborCells, + FullStencil, UniqueCellsView, ParticleRanges, + CellLookup, DenseLookup, UseDense, CellOccupied) + TargetLen = length(UniqueCellsView) + resize!(NeighborCellOffsets, TargetLen) + resize!(NeighborCellCounts, TargetLen) - if TargetLen > OriginalLen - @inbounds for Index in (OriginalLen + 1):TargetLen - NeighborCellLists[Index] = Int[] - end + MaxNeighbors = length(FullStencil) - 1 + RequiredLen = TargetLen * MaxNeighbors + if length(NeighborCells) < RequiredLen + resize!(NeighborCells, RequiredLen) end @inbounds for CellIndex in eachindex(UniqueCellsView) - Neighbors = NeighborCellLists[CellIndex] - empty!(Neighbors) - sizehint!(Neighbors, length(FullStencil) - 1) Cell = UniqueCellsView[CellIndex] + write_index = (CellIndex - 1) * MaxNeighbors + 1 + count = 0 for Offset in FullStencil NeighborCell = Cell + Offset - NeighborIndex = get(CellDict, NeighborCell, 0) - if NeighborIndex != 0 && NeighborIndex != CellIndex - StartIndex = ParticleRanges[NeighborIndex] - EndIndex = ParticleRanges[NeighborIndex + 1] - 1 - if StartIndex <= EndIndex - push!(Neighbors, NeighborIndex) - end + NeighborIndex = GetCellIndexDispatch(UseDense, DenseLookup, CellLookup, NeighborCell, 0) + if NeighborIndex != 0 && NeighborIndex != CellIndex && CellOccupied[NeighborIndex] + count += 1 + NeighborCells[write_index + count - 1] = NeighborIndex end end + NeighborCellOffsets[CellIndex] = write_index + NeighborCellCounts[CellIndex] = count end return nothing @@ -73,24 +244,30 @@ This builds a per-cell particle ordering buffer so cell ranges can be iterated without reordering the particle arrays. """ function UpdateNeighbors!(Particles, InverseCutOff, ParticleRanges, - UniqueCells, CellDict, ParticleOrder, CellOffsets) + UniqueCells, CellLookup, ParticleOrder, CellOffsets, + CellIndices, DenseLookup) ExtractCells!(Particles, InverseCutOff) Cells = @views Particles.Cells ParticleRanges[1] = 1 IndexCounter = 1 - empty!(CellDict) + ResetCellLookup!(CellLookup, length(Cells) * 8) fill!(CellOffsets, zero(eltype(CellOffsets))) + min_cell = Cells[1].I + max_cell = Cells[1].I @inbounds for Index in eachindex(Cells) Cell = Cells[Index] - CellIndex = get(CellDict, Cell, 0) - if CellIndex == 0 - IndexCounter += 1 - CellIndex = IndexCounter - CellDict[Cell] = CellIndex + for i in 1:length(Cell.I) + min_cell = Base.setindex(min_cell, min(min_cell[i], Cell.I[i]), i) + max_cell = Base.setindex(max_cell, max(max_cell[i], Cell.I[i]), i) + end + CellIndex, IsNew = GetOrInsertCellIndex!(CellLookup, Cell, IndexCounter + 1) + if IsNew + IndexCounter = CellIndex UniqueCells[CellIndex] = Cell end + CellIndices[Index] = CellIndex CellOffsets[CellIndex] += 1 end @@ -104,29 +281,63 @@ function UpdateNeighbors!(Particles, InverseCutOff, ParticleRanges, ParticleRanges[IndexCounter + 1] = RunningIndex @inbounds for Index in eachindex(Cells) - CellIndex = CellDict[Cells[Index]] + CellIndex = CellIndices[Index] TargetIndex = CellOffsets[CellIndex] ParticleOrder[TargetIndex] = Index CellOffsets[CellIndex] = TargetIndex + 1 end - return IndexCounter + dims = max_cell .- min_cell .+ 1 + dense_size = prod(dims) + use_dense = dense_size <= 8 * IndexCounter && dense_size <= 1_000_000 + if use_dense + DenseLookup.MinCell = SVector{length(min_cell), Int}(min_cell) + DenseLookup.Dims = SVector{length(dims), Int}(dims) + strides = ntuple(i -> i == 1 ? 1 : prod(dims[1:(i - 1)]), length(dims)) + DenseLookup.Strides = SVector{length(dims), Int}(strides) + if length(DenseLookup.IndexMap) < dense_size + resize!(DenseLookup.IndexMap, dense_size) + end + fill!(DenseLookup.IndexMap, 0) + @inbounds for CellIndex in 2:IndexCounter + Cell = UniqueCells[CellIndex] + linear = DenseIndex(DenseLookup.Dims, DenseLookup.Strides, DenseLookup.MinCell, Cell.I) + DenseLookup.IndexMap[linear] = CellIndex + end + end + + return IndexCounter, use_dense end function ComputeCellParticleCounts(ParticleRanges, CellCount) Counts = Vector{Int}(undef, CellCount) + ComputeCellParticleCounts!(Counts, ParticleRanges, CellCount) + return Counts +end + +function ComputeCellParticleCounts!(Counts, ParticleRanges, CellCount) + resize!(Counts, CellCount) @inbounds for Index in 1:CellCount Counts[Index] = ParticleRanges[Index + 1] - ParticleRanges[Index] end return Counts end -function ComputeCellNeighborCounts(ParticleRanges, NeighborCellLists, CellCount) +function ComputeCellNeighborCounts(ParticleRanges, NeighborCellOffsets, NeighborCellCounts, NeighborCells, CellCount) Counts = ComputeCellParticleCounts(ParticleRanges, CellCount) Neighbors = Vector{Int}(undef, CellCount) + ComputeCellNeighborCounts!(Neighbors, Counts, NeighborCellOffsets, NeighborCellCounts, NeighborCells, CellCount) + return Neighbors +end + +function ComputeCellNeighborCounts!(Neighbors, Counts, NeighborCellOffsets, NeighborCellCounts, NeighborCells, CellCount) + resize!(Neighbors, CellCount) @inbounds for Index in 1:CellCount NeighborTotal = 0 - for NeighborIndex in NeighborCellLists[Index] + count = NeighborCellCounts[Index] + start = NeighborCellOffsets[Index] + for j in 0:(count - 1) + NeighborIndex = NeighborCells[start + j] NeighborTotal += Counts[NeighborIndex] end Neighbors[Index] = max(Counts[Index] - 1, 0) + NeighborTotal @@ -134,6 +345,18 @@ function ComputeCellNeighborCounts(ParticleRanges, NeighborCellLists, CellCount) return Neighbors end +function UpdateCellCounts!(Scratch::NeighborListScratch, ParticleRanges, CellCount) + ResizeNeighborListScratch!(Scratch, CellCount) + Counts = Scratch.CellCounts + Occupied = Scratch.CellOccupied + @inbounds for Index in 1:CellCount + Count = ParticleRanges[Index + 1] - ParticleRanges[Index] + Counts[Index] = Count + Occupied[Index] = Count > 0 + end + return nothing +end + """ UpdateΔx!(Δx, posₙ⁺, pos)