diff --git a/src/SPHCellList.jl b/src/SPHCellList.jl index 3439367b..19943f6b 100644 --- a/src/SPHCellList.jl +++ b/src/SPHCellList.jl @@ -15,7 +15,7 @@ using ..OpenExternalPrograms using ..SPHKernels using ..SPHViscosityModels using ..SPHDensityDiffusionModels -using ..SPHNeighborList: BuildNeighborCellLists!, ComputeCellNeighborCounts, ComputeCellParticleCounts, ConstructStencil, ExtractCells!, FindCellIndex, MapFloor, UpdateNeighbors!, UpdateΔx! +using ..SPHNeighborList: BuildNeighborCellLists!, ComputeCellNeighborCounts, ComputeCellParticleCounts, ConstructStencil, ExtractCells!, MapFloor, SelectNeighborCellIndexType, UpdateNeighbors!, UpdateΔx! using Base.Threads: @threads using Bumper: @alloc, @no_escape @@ -29,7 +29,7 @@ using TimerOutputs: @timeit function NeighborLoopPerParticle!(SimDensityDiffusion::SDD, SimViscosity::SV, SimKernel, SimMetaData::SimulationMetaData{D,T,NoShifting,NoKernelOutput,B,L}, SimConstants, SimParticles, ParticleRanges, - CellListIndices, NeighborCellLists, dρdtI, + CellListIndices, NeighborCellListOffsets, NeighborCellListIndices, dρdtI, Acceleration, ∇Cᵢ, ∇◌rᵢ, AccelerationMax; Position = SimParticles.Position, @@ -46,7 +46,8 @@ using TimerOutputs: @timeit CellListIndex = CellListIndices[i] SameCellStart = ParticleRanges[CellListIndex] SameCellEnd = ParticleRanges[CellListIndex + 1] - 1 - NeighborCellIndices = NeighborCellLists[CellListIndex] + NeighborCellStart = Int(NeighborCellListOffsets[CellListIndex]) + NeighborCellEnd = Int(NeighborCellListOffsets[CellListIndex + 1]) - 1 @inbounds for j in SameCellStart:(i - 1) dρdt_acc, acc_acc = ComputeInteractionsPerParticle!( @@ -62,7 +63,8 @@ using TimerOutputs: @timeit Velocity, ParticleType, dρdt_acc, acc_acc, i, j, ) end - for NeighborIdx in NeighborCellIndices + for NeighborCellCursor in NeighborCellStart:NeighborCellEnd + NeighborIdx = Int(NeighborCellListIndices[NeighborCellCursor]) StartIndex_ = ParticleRanges[NeighborIdx] EndIndex_ = ParticleRanges[NeighborIdx + 1] - 1 @inbounds for j in StartIndex_:EndIndex_ @@ -85,7 +87,7 @@ using TimerOutputs: @timeit function NeighborLoopPerParticle!(SimDensityDiffusion::SDD, SimViscosity::SV, SimKernel, SimMetaData::SimulationMetaData{D,T,NoShifting,K,B,L}, SimConstants, SimParticles, ParticleRanges, - CellListIndices, NeighborCellLists, dρdtI, + CellListIndices, NeighborCellListOffsets, NeighborCellListIndices, dρdtI, Acceleration, ∇Cᵢ, ∇◌rᵢ, AccelerationMax; Position = SimParticles.Position, @@ -106,7 +108,8 @@ using TimerOutputs: @timeit CellListIndex = CellListIndices[i] SameCellStart = ParticleRanges[CellListIndex] SameCellEnd = ParticleRanges[CellListIndex + 1] - 1 - NeighborCellIndices = NeighborCellLists[CellListIndex] + NeighborCellStart = Int(NeighborCellListOffsets[CellListIndex]) + NeighborCellEnd = Int(NeighborCellListOffsets[CellListIndex + 1]) - 1 @inbounds for j in SameCellStart:(i - 1) dρdt_acc, acc_acc, kernel_acc, kernel_grad_acc = @@ -126,7 +129,8 @@ using TimerOutputs: @timeit kernel_grad_acc, i, j, ) end - for NeighborIdx in NeighborCellIndices + for NeighborCellCursor in NeighborCellStart:NeighborCellEnd + NeighborIdx = Int(NeighborCellListIndices[NeighborCellCursor]) StartIndex_ = ParticleRanges[NeighborIdx] EndIndex_ = ParticleRanges[NeighborIdx + 1] - 1 @inbounds for j in StartIndex_:EndIndex_ @@ -153,7 +157,7 @@ using TimerOutputs: @timeit function NeighborLoopPerParticle!(SimDensityDiffusion::SDD, SimViscosity::SV, SimKernel, SimMetaData::SimulationMetaData{D,T,S,NoKernelOutput,B,L}, SimConstants, SimParticles, ParticleRanges, - CellListIndices, NeighborCellLists, dρdtI, + CellListIndices, NeighborCellListOffsets, NeighborCellListIndices, dρdtI, Acceleration, ∇Cᵢ, ∇◌rᵢ, AccelerationMax; Position = SimParticles.Position, @@ -172,7 +176,8 @@ using TimerOutputs: @timeit CellListIndex = CellListIndices[i] SameCellStart = ParticleRanges[CellListIndex] SameCellEnd = ParticleRanges[CellListIndex + 1] - 1 - NeighborCellIndices = NeighborCellLists[CellListIndex] + NeighborCellStart = Int(NeighborCellListOffsets[CellListIndex]) + NeighborCellEnd = Int(NeighborCellListOffsets[CellListIndex + 1]) - 1 @inbounds for j in SameCellStart:(i - 1) dρdt_acc, acc_acc, shift_c_acc, shift_r_acc = @@ -192,7 +197,8 @@ using TimerOutputs: @timeit shift_r_acc, i, j, ) end - for NeighborIdx in NeighborCellIndices + for NeighborCellCursor in NeighborCellStart:NeighborCellEnd + NeighborIdx = Int(NeighborCellListIndices[NeighborCellCursor]) StartIndex_ = ParticleRanges[NeighborIdx] EndIndex_ = ParticleRanges[NeighborIdx + 1] - 1 @inbounds for j in StartIndex_:EndIndex_ @@ -219,7 +225,7 @@ using TimerOutputs: @timeit function NeighborLoopPerParticle!(SimDensityDiffusion::SDD, SimViscosity::SV, SimKernel, SimMetaData::SimulationMetaData{D,T,S,K,B,L}, SimConstants, SimParticles, ParticleRanges, - CellListIndices, NeighborCellLists, dρdtI, + CellListIndices, NeighborCellListOffsets, NeighborCellListIndices, dρdtI, Acceleration, ∇Cᵢ, ∇◌rᵢ, AccelerationMax; Position = SimParticles.Position, @@ -243,7 +249,8 @@ using TimerOutputs: @timeit CellListIndex = CellListIndices[i] SameCellStart = ParticleRanges[CellListIndex] SameCellEnd = ParticleRanges[CellListIndex + 1] - 1 - NeighborCellIndices = NeighborCellLists[CellListIndex] + NeighborCellStart = Int(NeighborCellListOffsets[CellListIndex]) + NeighborCellEnd = Int(NeighborCellListOffsets[CellListIndex + 1]) - 1 @inbounds for j in SameCellStart:(i - 1) dρdt_acc, acc_acc, kernel_acc, kernel_grad_acc, shift_c_acc, @@ -263,7 +270,8 @@ using TimerOutputs: @timeit kernel_grad_acc, shift_c_acc, shift_r_acc, i, j, ) end - for NeighborIdx in NeighborCellIndices + for NeighborCellCursor in NeighborCellStart:NeighborCellEnd + NeighborIdx = Int(NeighborCellListIndices[NeighborCellCursor]) StartIndex_ = ParticleRanges[NeighborIdx] EndIndex_ = ParticleRanges[NeighborIdx + 1] - 1 @inbounds for j in StartIndex_:EndIndex_ @@ -289,17 +297,15 @@ using TimerOutputs: @timeit return nothing end - f(SimKernel, GhostPoint) = CartesianIndex(map(x -> MapFloor(x, SimKernel.H⁻¹), Tuple(GhostPoint))) + @inline GhostCellIndex(SimKernel, GhostPoint) = CartesianIndex(map(x -> MapFloor(x, SimKernel.H⁻¹), Tuple(GhostPoint))) function NeighborLoopMDBC!(SimKernel, SimMetaData::SimulationMetaData{Dimensions, FloatType, SMode, KMode, BMode, LMode}, - SimConstants, ParticleRanges, UniqueCellsView, - SimParticles, bᵧ, Aᵧ) where {Dimensions, FloatType, SMode, KMode, BMode, LMode} + SimConstants, FullStencil, ParticleRanges, + CellIndexMap::Dict{CellType, IndexType}, SimParticles, bᵧ, Aᵧ) where {Dimensions, FloatType, SMode, KMode, BMode, LMode, CellType, IndexType<:Integer} - @unpack Position, Density, GhostPoints, GhostNormals = SimParticles + @unpack Position, Density, GhostPoints = SimParticles ParticleType = SimParticles.Type - FullStencil = ConstructStencil(Val(Dimensions)) - @inbounds @threads for iter in eachindex(GhostPoints) GhostPoint = GhostPoints[iter] @@ -309,14 +315,13 @@ using TimerOutputs: @timeit A_acc = zero(Aᵧ[iter]) # an SMatrix{D+1,D+1,FloatType} # compute and accumulate into the locals - GhostCellIndex = f(SimKernel, GhostPoints[iter]) - @inbounds for offset ∈ FullStencil - SCellIndex = GhostCellIndex + offset - - # 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 = FindCellIndex(UniqueCellsView, SCellIndex) + CellIndex = GhostCellIndex(SimKernel, GhostPoint) + @inbounds for Offset ∈ FullStencil + NeighborIdxValue = get(CellIndexMap, CellIndex + Offset, zero(IndexType)) + if iszero(NeighborIdxValue) + continue + end + NeighborIdx = Int(NeighborIdxValue) StartIndex_ = ParticleRanges[NeighborIdx] EndIndex_ = ParticleRanges[NeighborIdx + 1] - 1 @@ -628,14 +633,13 @@ using TimerOutputs: @timeit function ApplyMDBCBeforeHalf!(SimMetaData::SimulationMetaData{D,T,S,K,SimpleMDBC,L}, SimKernel, SimConstants, SimParticles, - ParticleRanges, UniqueCells + FullStencil, ParticleRanges, CellIndexMap ) where {D,T,S<:ShiftingMode,K<:KernelOutputMode,L<:LogMode} @no_escape begin DimensionsPlus = D + 1 bᵧ = @alloc(SVector{DimensionsPlus, T}, length(SimParticles.Position)) Aᵧ = @alloc(SMatrix{DimensionsPlus, DimensionsPlus, T, DimensionsPlus*DimensionsPlus}, length(SimParticles.Position)) - UniqueCellsView = view(UniqueCells, 1:SimMetaData.IndexCounter) - NeighborLoopMDBC!(SimKernel, SimMetaData, SimConstants, ParticleRanges, UniqueCellsView, SimParticles, bᵧ, Aᵧ) + NeighborLoopMDBC!(SimKernel, SimMetaData, SimConstants, FullStencil, ParticleRanges, CellIndexMap, SimParticles, bᵧ, Aᵧ) ApplyMDBCCorrection(SimConstants, SimParticles, bᵧ, Aᵧ) end @@ -696,7 +700,7 @@ using TimerOutputs: @timeit SimConstants, SimParticles, FullStencil, ParticleRanges, UniqueCells, CellListIndices, SortingScratchSpace, - NeighborCellLists, dρdtI, Velocityₙ⁺, + NeighborCellListOffsets, NeighborCellListIndices, NeighborCellIndexMap, dρdtI, Velocityₙ⁺, Positionₙ⁺, ρₙ⁺, ∇Cᵢ, ∇◌rᵢ, MotionDefinition::Union{ Nothing, @@ -728,15 +732,15 @@ using TimerOutputs: @timeit SimMetaData.IndexCounter = UpdateNeighbors!(SimParticles, SimKernel.H⁻¹, SortingScratchSpace, ParticleRanges, UniqueCells, CellListIndices) UniqueCellsView = view(UniqueCells, 1:SimMetaData.IndexCounter) - BuildNeighborCellLists!(NeighborCellLists, FullStencil, UniqueCellsView, ParticleRanges) + BuildNeighborCellLists!(NeighborCellListOffsets, NeighborCellListIndices, FullStencil, UniqueCellsView, ParticleRanges, NeighborCellIndexMap) if TimeSteppingMode isa SingleNeighborTimeStepping @timeit SimMetaData.HourGlass "00 Init Pressure" Pressure!(SimParticles.Pressure, SimParticles.Density, SimConstants) - @timeit SimMetaData.HourGlass "00a Init MDBC" ApplyMDBCBeforeHalf!(SimMetaData, SimKernel, SimConstants, SimParticles, ParticleRanges, UniqueCells) + @timeit SimMetaData.HourGlass "00a Init MDBC" ApplyMDBCBeforeHalf!(SimMetaData, SimKernel, SimConstants, SimParticles, FullStencil, ParticleRanges, NeighborCellIndexMap) @timeit SimMetaData.HourGlass "00b Init NeighborLoop" NeighborLoopPerParticle!( SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData, SimConstants, SimParticles, ParticleRanges, CellListIndices, - NeighborCellLists, dρdtI, SimParticles.Acceleration, ∇Cᵢ, ∇◌rᵢ, AccelerationMax, + NeighborCellListOffsets, NeighborCellListIndices, dρdtI, SimParticles.Acceleration, ∇Cᵢ, ∇◌rᵢ, AccelerationMax, ) end @@ -759,7 +763,7 @@ using TimerOutputs: @timeit @timeit SimMetaData.HourGlass "01a Actual Calculate IndexCounter" SimMetaData.IndexCounter = UpdateNeighbors!(SimParticles, SimKernel.H⁻¹, SortingScratchSpace, ParticleRanges, UniqueCells, CellListIndices) SimMetaData.Δx = zero(eltype(dρdtI)) UniqueCellsView = view(UniqueCells, 1:SimMetaData.IndexCounter) - BuildNeighborCellLists!(NeighborCellLists, FullStencil, UniqueCellsView, ParticleRanges) + BuildNeighborCellLists!(NeighborCellListOffsets, NeighborCellListIndices, FullStencil, UniqueCellsView, ParticleRanges, NeighborCellIndexMap) end end @@ -767,12 +771,12 @@ using TimerOutputs: @timeit if TimeSteppingMode isa SymplecticTimeStepping @timeit SimMetaData.HourGlass "02 Pressure" Pressure!(SimParticles.Pressure, SimParticles.Density, SimConstants) - @timeit SimMetaData.HourGlass "03 Apply MDBC before Half TimeStep" ApplyMDBCBeforeHalf!(SimMetaData, SimKernel, SimConstants, SimParticles, ParticleRanges, UniqueCells) + @timeit SimMetaData.HourGlass "03 Apply MDBC before Half TimeStep" ApplyMDBCBeforeHalf!(SimMetaData, SimKernel, SimConstants, SimParticles, FullStencil, ParticleRanges, NeighborCellIndexMap) @timeit SimMetaData.HourGlass "04 First NeighborLoop" NeighborLoopPerParticle!( SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData, SimConstants, SimParticles, ParticleRanges, CellListIndices, - NeighborCellLists, dρdtI, SimParticles.Acceleration, ∇Cᵢ, ∇◌rᵢ, AccelerationMax, + NeighborCellListOffsets, NeighborCellListIndices, dρdtI, SimParticles.Acceleration, ∇Cᵢ, ∇◌rᵢ, AccelerationMax, ) @timeit SimMetaData.HourGlass "05 Update To Half TimeStep" HalfTimeStep(SimMetaData, SimConstants, SimParticles, Positionₙ⁺, Velocityₙ⁺, ρₙ⁺, dρdtI, dt₂) @@ -785,13 +789,13 @@ using TimerOutputs: @timeit @timeit SimMetaData.HourGlass "08 Second NeighborLoop" NeighborLoopPerParticle!( SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData, SimConstants, SimParticles, ParticleRanges, CellListIndices, - NeighborCellLists, dρdtI, SimParticles.Acceleration, ∇Cᵢ, ∇◌rᵢ, AccelerationMax, + NeighborCellListOffsets, NeighborCellListIndices, dρdtI, SimParticles.Acceleration, ∇Cᵢ, ∇◌rᵢ, AccelerationMax, Position = Positionₙ⁺, Density = ρₙ⁺, Velocity = Velocityₙ⁺, ) else - @timeit SimMetaData.HourGlass "02 Apply MDBC before Half TimeStep" ApplyMDBCBeforeHalf!(SimMetaData, SimKernel, SimConstants, SimParticles, ParticleRanges, UniqueCells) + @timeit SimMetaData.HourGlass "02 Apply MDBC before Half TimeStep" ApplyMDBCBeforeHalf!(SimMetaData, SimKernel, SimConstants, SimParticles, FullStencil, ParticleRanges, NeighborCellIndexMap) @timeit SimMetaData.HourGlass "03 Update To Half TimeStep" HalfTimeStep(SimMetaData, SimConstants, SimParticles, Positionₙ⁺, Velocityₙ⁺, ρₙ⁺, dρdtI, dt₂) @@ -803,7 +807,7 @@ using TimerOutputs: @timeit @timeit SimMetaData.HourGlass "06 NeighborLoop" NeighborLoopPerParticle!( SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData, SimConstants, SimParticles, ParticleRanges, CellListIndices, - NeighborCellLists, dρdtI, SimParticles.Acceleration, ∇Cᵢ, ∇◌rᵢ, AccelerationMax, + NeighborCellListOffsets, NeighborCellListIndices, dρdtI, SimParticles.Acceleration, ∇Cᵢ, ∇◌rᵢ, AccelerationMax, Position = Positionₙ⁺, Density = ρₙ⁺, Velocity = Velocityₙ⁺, @@ -911,7 +915,10 @@ using TimerOutputs: @timeit UniqueCells = zeros(CartesianIndex{Dimensions}, NumberOfPoints) CellListIndices = zeros(Int, NumberOfPoints) FullStencil = ConstructStencil(Val(Dimensions)) - NeighborCellLists = [Int[] for _ in 1:length(UniqueCells)] + NeighborIndexType = SelectNeighborCellIndexType(length(UniqueCells), FullStencil) + NeighborCellListOffsets = zeros(NeighborIndexType, length(UniqueCells) + 1) + NeighborCellListIndices = NeighborIndexType[] + NeighborCellIndexMap = Dict{CartesianIndex{Dimensions}, NeighborIndexType}() _, SortingScratchSpace = Base.Sort.make_scratch(nothing, eltype(SimParticles), NumberOfPoints) output = SetupVTKOutput(SimMetaData, SimParticles, SimKernel, Dimensions) @@ -930,7 +937,8 @@ using TimerOutputs: @timeit ) cell_neighbor_counts = ComputeCellNeighborCounts( ParticleRanges, - NeighborCellLists, + NeighborCellListOffsets, + NeighborCellListIndices, SimMetaData.IndexCounter, ) end @@ -952,7 +960,7 @@ using TimerOutputs: @timeit SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData, SimConstants, SimParticles, FullStencil, ParticleRanges, UniqueCells, CellListIndices, SortingScratchSpace, - NeighborCellLists, dρdtI, Velocityₙ⁺, Positionₙ⁺, ρₙ⁺, + NeighborCellListOffsets, NeighborCellListIndices, NeighborCellIndexMap, dρdtI, Velocityₙ⁺, Positionₙ⁺, ρₙ⁺, ∇Cᵢ, ∇◌rᵢ, MotionDefinition, ) push!(SimMetaData.TimeSteps, SimMetaData.CurrentTimeStep) @@ -971,7 +979,8 @@ using TimerOutputs: @timeit ) cell_neighbor_counts = ComputeCellNeighborCounts( ParticleRanges, - NeighborCellLists, + NeighborCellListOffsets, + NeighborCellListIndices, length(UniqueCellsView), ) end diff --git a/src/SPHNeighborList.jl b/src/SPHNeighborList.jl index cab928a8..be046cf3 100644 --- a/src/SPHNeighborList.jl +++ b/src/SPHNeighborList.jl @@ -4,6 +4,19 @@ export ConstructStencil, ExtractCells!, UpdateNeighbors!, BuildNeighborCellLists using StaticArrays +const NeighborCellIndexType = UInt32 + +function SelectNeighborCellIndexType(CellCapacity, FullStencil) + MaxNeighborCount = max(length(FullStencil) - 1, 0) + MaxStoredIndex = max(CellCapacity, CellCapacity * MaxNeighborCount + 1) + if MaxStoredIndex <= typemax(UInt16) + return UInt16 + elseif MaxStoredIndex <= typemax(UInt32) + return UInt32 + end + return UInt64 +end + function ConstructStencil(V::Val{d}) where d return CartesianIndices(ntuple(_ -> -1:1, V)) end @@ -21,26 +34,49 @@ end end function BuildNeighborCellLists!(NeighborCellLists, FullStencil, UniqueCellsView, ParticleRanges) + CellIndexMap = Dict{eltype(UniqueCellsView), Int}() + return BuildNeighborCellLists!(NeighborCellLists, FullStencil, UniqueCellsView, ParticleRanges, CellIndexMap) +end + +@inline function CompressIndex(::Type{IndexType}, Index) where {IndexType<:Integer} + @assert Index <= typemax(IndexType) + return IndexType(Index) +end + +function BuildNeighborCellLists!(NeighborCellLists, FullStencil, UniqueCellsView, ParticleRanges, + CellIndexMap::Dict{CellType, IndexType}) where {CellType, IndexType<:Integer} TargetLen = length(UniqueCellsView) OriginalLen = length(NeighborCellLists) resize!(NeighborCellLists, TargetLen) + MaxNeighborCount = max(length(FullStencil) - 1, 0) if TargetLen > OriginalLen @inbounds for Index in (OriginalLen + 1):TargetLen NeighborCellLists[Index] = Int[] + sizehint!(NeighborCellLists[Index], MaxNeighborCount) + end + end + + empty!(CellIndexMap) + sizehint!(CellIndexMap, TargetLen) + @inbounds for CellIndex in eachindex(UniqueCellsView) + if ParticleRanges[CellIndex] < ParticleRanges[CellIndex + 1] + CellIndexMap[UniqueCellsView[CellIndex]] = CompressIndex(IndexType, CellIndex) end end @inbounds for CellIndex in eachindex(UniqueCellsView) Neighbors = NeighborCellLists[CellIndex] empty!(Neighbors) + if ParticleRanges[CellIndex] >= ParticleRanges[CellIndex + 1] + continue + end Cell = UniqueCellsView[CellIndex] for Offset in FullStencil NeighborCell = Cell + Offset - NeighborIndex = FindCellIndex(UniqueCellsView, NeighborCell) - StartIndex = ParticleRanges[NeighborIndex] - EndIndex = ParticleRanges[NeighborIndex + 1] - 1 - if StartIndex <= EndIndex && NeighborIndex != CellIndex + NeighborIndexValue = get(CellIndexMap, NeighborCell, zero(IndexType)) + NeighborIndex = Int(NeighborIndexValue) + if !iszero(NeighborIndex) && NeighborIndex != CellIndex push!(Neighbors, NeighborIndex) end end @@ -49,6 +85,42 @@ function BuildNeighborCellLists!(NeighborCellLists, FullStencil, UniqueCellsView return nothing end +function BuildNeighborCellLists!(NeighborCellListOffsets, NeighborCellListIndices, + FullStencil, UniqueCellsView, ParticleRanges, + CellIndexMap::Dict{CellType, IndexType}) where {CellType, IndexType<:Integer} + TargetLen = length(UniqueCellsView) + resize!(NeighborCellListOffsets, TargetLen + 1) + MaxNeighborCount = max(length(FullStencil) - 1, 0) + + empty!(CellIndexMap) + sizehint!(CellIndexMap, TargetLen) + @inbounds for CellIndex in eachindex(UniqueCellsView) + if ParticleRanges[CellIndex] < ParticleRanges[CellIndex + 1] + CellIndexMap[UniqueCellsView[CellIndex]] = CompressIndex(IndexType, CellIndex) + end + end + + empty!(NeighborCellListIndices) + sizehint!(NeighborCellListIndices, TargetLen * MaxNeighborCount) + NeighborCellListOffsets[1] = CompressIndex(eltype(NeighborCellListOffsets), 1) + @inbounds for CellIndex in eachindex(UniqueCellsView) + NeighborCellListOffsets[CellIndex] = CompressIndex(eltype(NeighborCellListOffsets), length(NeighborCellListIndices) + 1) + if ParticleRanges[CellIndex] < ParticleRanges[CellIndex + 1] + Cell = UniqueCellsView[CellIndex] + for Offset in FullStencil + NeighborIndexValue = get(CellIndexMap, Cell + Offset, zero(IndexType)) + NeighborIndex = Int(NeighborIndexValue) + if !iszero(NeighborIndex) && NeighborIndex != CellIndex + push!(NeighborCellListIndices, CompressIndex(eltype(NeighborCellListIndices), NeighborIndex)) + end + end + end + end + NeighborCellListOffsets[TargetLen + 1] = CompressIndex(eltype(NeighborCellListOffsets), length(NeighborCellListIndices) + 1) + + return nothing +end + """ Extracts the cells for each particle based on their positions and the inverse cutoff value. @@ -135,6 +207,19 @@ function ComputeCellNeighborCounts(ParticleRanges, NeighborCellLists, CellCount) return Neighbors end +function ComputeCellNeighborCounts(ParticleRanges, NeighborCellListOffsets, NeighborCellListIndices, CellCount) + Counts = ComputeCellParticleCounts(ParticleRanges, CellCount) + Neighbors = Vector{Int}(undef, CellCount) + @inbounds for Index in 1:CellCount + NeighborTotal = 0 + for Cursor in Int(NeighborCellListOffsets[Index]):(Int(NeighborCellListOffsets[Index + 1]) - 1) + NeighborTotal += Counts[Int(NeighborCellListIndices[Cursor])] + end + Neighbors[Index] = max(Counts[Index] - 1, 0) + NeighborTotal + end + return Neighbors +end + """ UpdateΔx!(Δx, posₙ⁺, pos)