From b51c1c507b6f92623ef5fe58149dd94c6a86a871 Mon Sep 17 00:00:00 2001 From: AhmedSalih3d <36305327+AhmedSalih3d@users.noreply.github.com> Date: Thu, 15 Jan 2026 20:55:51 +0100 Subject: [PATCH 1/6] Replace cell dict with lookup arrays --- src/SPHCellList.jl | 46 ++++++++++---------- src/SPHNeighborList.jl | 98 ++++++++++++++++++++++++++++++++++++++---- 2 files changed, 113 insertions(+), 31 deletions(-) diff --git a/src/SPHCellList.jl b/src/SPHCellList.jl index 19ddae69..89111a3d 100644 --- a/src/SPHCellList.jl +++ b/src/SPHCellList.jl @@ -18,7 +18,7 @@ using ..OpenExternalPrograms using ..SPHKernels using ..SPHViscosityModels using ..SPHDensityDiffusionModels -using ..SPHNeighborList: BuildNeighborCellLists!, ComputeCellNeighborCounts, ComputeCellParticleCounts, ConstructStencil, ExtractCells!, MapFloor, UpdateNeighbors!, UpdateΔx! +using ..SPHNeighborList: BuildNeighborCellLists!, CellLookup, ComputeCellNeighborCounts, ComputeCellParticleCounts, ConstructStencil, ExtractCells!, GetCellIndex, InitializeCellLookup, MapFloor, UpdateNeighbors!, UpdateΔx! using StaticArrays using StructArrays: StructArray, foreachfield @@ -35,7 +35,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, + CellLookup, NeighborCellLists, dρdtI, Acceleration, ∇Cᵢ, ∇◌rᵢ, max_visc = nothing, min_dt_force = nothing; @@ -52,7 +52,7 @@ using LinearAlgebra dρdt_acc = zero(dρdtI[i]) acc_acc = zero(Acceleration[i]) CellIndex = Cells[i] - CellListIndex = get(CellDict, CellIndex, 1) + CellListIndex = GetCellIndex(CellLookup, CellIndex, 1) SameCellStart = ParticleRanges[CellListIndex] SameCellEnd = ParticleRanges[CellListIndex + 1] - 1 NeighborCellIndices = NeighborCellLists[CellListIndex] @@ -91,7 +91,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, + CellLookup, NeighborCellLists, dρdtI, Acceleration, ∇Cᵢ, ∇◌rᵢ, max_visc = nothing, min_dt_force = nothing; @@ -111,7 +111,7 @@ using LinearAlgebra kernel_acc = zero(Kernel[i]) kernel_grad_acc = zero(KernelGradient[i]) CellIndex = Cells[i] - CellListIndex = get(CellDict, CellIndex, 1) + CellListIndex = GetCellIndex(CellLookup, CellIndex, 1) SameCellStart = ParticleRanges[CellListIndex] SameCellEnd = ParticleRanges[CellListIndex + 1] - 1 NeighborCellIndices = NeighborCellLists[CellListIndex] @@ -156,7 +156,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, + CellLookup, NeighborCellLists, dρdtI, Acceleration, ∇Cᵢ, ∇◌rᵢ, max_visc = nothing, min_dt_force = nothing; @@ -175,7 +175,7 @@ using LinearAlgebra shift_c_acc = zero(∇Cᵢ[i]) shift_r_acc = zero(∇◌rᵢ[i]) CellIndex = Cells[i] - CellListIndex = get(CellDict, CellIndex, 1) + CellListIndex = GetCellIndex(CellLookup, CellIndex, 1) SameCellStart = ParticleRanges[CellListIndex] SameCellEnd = ParticleRanges[CellListIndex + 1] - 1 NeighborCellIndices = NeighborCellLists[CellListIndex] @@ -220,7 +220,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, + CellLookup, NeighborCellLists, dρdtI, Acceleration, ∇Cᵢ, ∇◌rᵢ, max_visc = nothing, min_dt_force = nothing; @@ -243,7 +243,7 @@ using LinearAlgebra shift_c_acc = zero(∇Cᵢ[i]) shift_r_acc = zero(∇◌rᵢ[i]) CellIndex = Cells[i] - CellListIndex = get(CellDict, CellIndex, 1) + CellListIndex = GetCellIndex(CellLookup, CellIndex, 1) SameCellStart = ParticleRanges[CellListIndex] SameCellEnd = ParticleRanges[CellListIndex + 1] - 1 NeighborCellIndices = NeighborCellLists[CellListIndex] @@ -291,7 +291,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, Position, Density, GhostPoints, GhostNormals, ParticleType, bᵧ, Aᵧ; ParticleOrder) where {Dimensions, FloatType, SMode, KMode, BMode, LMode} @@ -314,7 +314,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 = GetCellIndex(CellLookup, SCellIndex, 1) StartIndex_ = ParticleRanges[NeighborIdx] EndIndex_ = ParticleRanges[NeighborIdx + 1] - 1 @@ -632,7 +632,7 @@ using LinearAlgebra function ApplyMDBCBeforeHalf!(SimMetaData::SimulationMetaData{D,T,S,K,SimpleMDBC,L}, SimKernel, SimConstants, SimParticles, - ParticleRanges, CellDict, Position, Density, + ParticleRanges, CellLookup, Position, Density, GhostPoints, GhostNormals, ParticleType; ParticleOrder ) where {D,T,S<:ShiftingMode,K<:KernelOutputMode,L<:LogMode} @@ -640,7 +640,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, Position, Density, GhostPoints, GhostNormals, ParticleType, bᵧ, Aᵧ; ParticleOrder = ParticleOrder) ApplyMDBCCorrection(SimConstants, SimParticles, bᵧ, Aᵧ) end @@ -697,7 +697,7 @@ using LinearAlgebra @inbounds function SimulationLoop(SimDensityDiffusion::SDD, SimViscosity::SV, SimKernel, SimMetaData::SimulationMetaData{Dimensions, FloatType, SMode, KMode, BMode, LMode}, SimConstants, SimParticles, FullStencil, - ParticleRanges, UniqueCells, CellDict, + ParticleRanges, UniqueCells, CellLookup, ParticleOrder, CellOffsets, NeighborCellLists, dρdtI, Velocityₙ⁺, Positionₙ⁺, ρₙ⁺, ∇Cᵢ, ∇◌rᵢ, @@ -747,10 +747,10 @@ 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" SimMetaData.IndexCounter = UpdateNeighbors!(SimParticles, SimKernel.H⁻¹, ParticleRanges, UniqueCells, CellLookup, ParticleOrder, CellOffsets) SimMetaData.Δx = zero(eltype(dρdtI)) UniqueCellsView = view(UniqueCells, 1:SimMetaData.IndexCounter) - BuildNeighborCellLists!(NeighborCellLists, FullStencil, UniqueCellsView, ParticleRanges, CellDict) + BuildNeighborCellLists!(NeighborCellLists, FullStencil, UniqueCellsView, ParticleRanges, CellLookup) end end @@ -760,20 +760,20 @@ 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, 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, + SimConstants, SimParticles, ParticleRanges, CellLookup, NeighborCellLists, dρdtI, Acceleration, ∇Cᵢ, ∇◌rᵢ, ParticleOrder = ParticleOrder, ) else @timeit SimMetaData.HourGlass "04 First NeighborLoop" NeighborLoopPerParticle!( SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData, - SimConstants, SimParticles, ParticleRanges, CellDict, + SimConstants, SimParticles, ParticleRanges, CellLookup, NeighborCellLists, dρdtI, Acceleration, ∇Cᵢ, ∇◌rᵢ, ParticleOrder = ParticleOrder, ) @@ -791,7 +791,7 @@ 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, + SimConstants, SimParticles, ParticleRanges, CellLookup, NeighborCellLists, dρdtI, Acceleration, ∇Cᵢ, ∇◌rᵢ, max_visc, min_dt_force, ParticleOrder = ParticleOrder, @@ -802,7 +802,7 @@ using LinearAlgebra else @timeit SimMetaData.HourGlass "08 Second NeighborLoop" NeighborLoopPerParticle!( SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData, - SimConstants, SimParticles, ParticleRanges, CellDict, + SimConstants, SimParticles, ParticleRanges, CellLookup, NeighborCellLists, dρdtI, Acceleration, ∇Cᵢ, ∇◌rᵢ, max_visc, min_dt_force, ParticleOrder = ParticleOrder, @@ -867,7 +867,7 @@ 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) FullStencil = ConstructStencil(Val(Dimensions)) NeighborCellLists = [Int[] for _ in 1:length(UniqueCells)] ParticleOrder = zeros(Int, NumberOfPoints) @@ -909,7 +909,7 @@ using LinearAlgebra @timeit SimMetaData.HourGlass "00 SimulationLoop" SimulationLoop( SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData, SimConstants, SimParticles, FullStencil, ParticleRanges, - UniqueCells, CellDict, ParticleOrder, CellOffsets, + UniqueCells, CellLookup, ParticleOrder, CellOffsets, NeighborCellLists, dρdtI, Velocityₙ⁺, Positionₙ⁺, ρₙ⁺, ∇Cᵢ, ∇◌rᵢ, MotionDefinition, ) diff --git a/src/SPHNeighborList.jl b/src/SPHNeighborList.jl index b94fc07e..5ba90da1 100644 --- a/src/SPHNeighborList.jl +++ b/src/SPHNeighborList.jl @@ -1,15 +1,97 @@ module SPHNeighborList export ConstructStencil, ExtractCells!, UpdateNeighbors!, - BuildNeighborCellLists!, ComputeCellParticleCounts, ComputeCellNeighborCounts, UpdateΔx! + BuildNeighborCellLists!, ComputeCellParticleCounts, ComputeCellNeighborCounts, UpdateΔx!, + CellLookup, InitializeCellLookup, ResetCellLookup!, GetCellIndex, SetCellIndex! using StaticArrays +mutable struct CellLookup{D} + Keys::Vector{CartesianIndex{D}} + Values::Vector{Int} + Filled::Vector{Bool} + Mask::Int +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), fill(false, size), size - 1) +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.Filled = fill(false, size) + Lookup.Mask = size - 1 + else + fill!(Lookup.Filled, false) + end + return nothing +end + +@inline function HashCellIndex(Cell::CartesianIndex{D}) where D + h = UInt(0x9e3779b97f4a7c15) + @inbounds for i in 1:D + v = UInt(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)) + 1) + @inbounds for _ in 0:mask + if !Lookup.Filled[index] + return default + elseif Lookup.Keys[index] == Cell + return Lookup.Values[index] + end + index = (index & mask) + 1 + 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)) + 1) + @inbounds for _ in 0:mask + if !Lookup.Filled[index] + Lookup.Filled[index] = true + Lookup.Keys[index] = Cell + Lookup.Values[index] = value + return nothing + elseif Lookup.Keys[index] == Cell + Lookup.Values[index] = value + return nothing + end + index = (index & mask) + 1 + end + return nothing +end + function ConstructStencil(V::Val{d}) where d return CartesianIndices(ntuple(_ -> -1:1, V)) end -function BuildNeighborCellLists!(NeighborCellLists, FullStencil, UniqueCellsView, ParticleRanges, CellDict) +function BuildNeighborCellLists!(NeighborCellLists, FullStencil, UniqueCellsView, ParticleRanges, CellLookup) TargetLen = length(UniqueCellsView) OriginalLen = length(NeighborCellLists) resize!(NeighborCellLists, TargetLen) @@ -27,7 +109,7 @@ function BuildNeighborCellLists!(NeighborCellLists, FullStencil, UniqueCellsView Cell = UniqueCellsView[CellIndex] for Offset in FullStencil NeighborCell = Cell + Offset - NeighborIndex = get(CellDict, NeighborCell, 0) + NeighborIndex = GetCellIndex(CellLookup, NeighborCell, 0) if NeighborIndex != 0 && NeighborIndex != CellIndex StartIndex = ParticleRanges[NeighborIndex] EndIndex = ParticleRanges[NeighborIndex + 1] - 1 @@ -73,22 +155,22 @@ 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) ExtractCells!(Particles, InverseCutOff) Cells = @views Particles.Cells ParticleRanges[1] = 1 IndexCounter = 1 - empty!(CellDict) + ResetCellLookup!(CellLookup, length(Cells) * 2) fill!(CellOffsets, zero(eltype(CellOffsets))) @inbounds for Index in eachindex(Cells) Cell = Cells[Index] - CellIndex = get(CellDict, Cell, 0) + CellIndex = GetCellIndex(CellLookup, Cell, 0) if CellIndex == 0 IndexCounter += 1 CellIndex = IndexCounter - CellDict[Cell] = CellIndex + SetCellIndex!(CellLookup, Cell, CellIndex) UniqueCells[CellIndex] = Cell end CellOffsets[CellIndex] += 1 @@ -104,7 +186,7 @@ function UpdateNeighbors!(Particles, InverseCutOff, ParticleRanges, ParticleRanges[IndexCounter + 1] = RunningIndex @inbounds for Index in eachindex(Cells) - CellIndex = CellDict[Cells[Index]] + CellIndex = GetCellIndex(CellLookup, Cells[Index], 0) TargetIndex = CellOffsets[CellIndex] ParticleOrder[TargetIndex] = Index CellOffsets[CellIndex] = TargetIndex + 1 From 16bfe3427c56e5a83183c6b25df07732c8b09734 Mon Sep 17 00:00:00 2001 From: AhmedSalih3d <36305327+AhmedSalih3d@users.noreply.github.com> Date: Thu, 15 Jan 2026 21:02:37 +0100 Subject: [PATCH 2/6] Fix cell lookup hashing for negatives --- src/SPHNeighborList.jl | 30 ++++++++++++++++-------------- 1 file changed, 16 insertions(+), 14 deletions(-) diff --git a/src/SPHNeighborList.jl b/src/SPHNeighborList.jl index 5ba90da1..9656c21a 100644 --- a/src/SPHNeighborList.jl +++ b/src/SPHNeighborList.jl @@ -49,7 +49,7 @@ end @inline function HashCellIndex(Cell::CartesianIndex{D}) where D h = UInt(0x9e3779b97f4a7c15) @inbounds for i in 1:D - v = UInt(Cell.I[i]) + v = unsigned(Cell.I[i]) h ⊻= v + 0x9e3779b97f4a7c15 + (h << 6) + (h >> 2) end return h @@ -57,32 +57,34 @@ end @inline function GetCellIndex(Lookup::CellLookup{D}, Cell::CartesianIndex{D}, default::Int) where D mask = Lookup.Mask - index = Int((HashCellIndex(Cell) & UInt(mask)) + 1) + index = Int(HashCellIndex(Cell) & UInt(mask)) @inbounds for _ in 0:mask - if !Lookup.Filled[index] + slot = index + 1 + if !Lookup.Filled[slot] return default - elseif Lookup.Keys[index] == Cell - return Lookup.Values[index] + elseif Lookup.Keys[slot] == Cell + return Lookup.Values[slot] end - index = (index & mask) + 1 + 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)) + 1) + index = Int(HashCellIndex(Cell) & UInt(mask)) @inbounds for _ in 0:mask - if !Lookup.Filled[index] - Lookup.Filled[index] = true - Lookup.Keys[index] = Cell - Lookup.Values[index] = value + slot = index + 1 + if !Lookup.Filled[slot] + Lookup.Filled[slot] = true + Lookup.Keys[slot] = Cell + Lookup.Values[slot] = value return nothing - elseif Lookup.Keys[index] == Cell - Lookup.Values[index] = value + elseif Lookup.Keys[slot] == Cell + Lookup.Values[slot] = value return nothing end - index = (index & mask) + 1 + index = (index + 1) & mask end return nothing end From 4cb4d9886f58d4f543a18a0dd209f7e1dd682a75 Mon Sep 17 00:00:00 2001 From: AhmedSalih3d <36305327+AhmedSalih3d@users.noreply.github.com> Date: Thu, 15 Jan 2026 21:10:05 +0100 Subject: [PATCH 3/6] Speed up cell lookup probing --- src/SPHNeighborList.jl | 50 +++++++++++++++++++++++++++++++----------- 1 file changed, 37 insertions(+), 13 deletions(-) diff --git a/src/SPHNeighborList.jl b/src/SPHNeighborList.jl index 9656c21a..8fbb9920 100644 --- a/src/SPHNeighborList.jl +++ b/src/SPHNeighborList.jl @@ -9,8 +9,9 @@ using StaticArrays mutable struct CellLookup{D} Keys::Vector{CartesianIndex{D}} Values::Vector{Int} - Filled::Vector{Bool} + Stamps::Vector{UInt32} Mask::Int + Epoch::UInt32 end @inline function ZeroCell(::Val{D}) where D @@ -28,7 +29,7 @@ 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), fill(false, size), size - 1) + return CellLookup{D}(fill(zero_cell, size), zeros(Int, size), zeros(UInt32, size), size - 1, UInt32(1)) end function ResetCellLookup!(Lookup::CellLookup{D}, capacity::Int) where D @@ -38,10 +39,15 @@ function ResetCellLookup!(Lookup::CellLookup{D}, capacity::Int) where D zero_cell = ZeroCell(Val(D)) Lookup.Keys = fill(zero_cell, size) Lookup.Values = zeros(Int, size) - Lookup.Filled = fill(false, size) + Lookup.Stamps = zeros(UInt32, size) Lookup.Mask = size - 1 + Lookup.Epoch = UInt32(1) else - fill!(Lookup.Filled, false) + Lookup.Epoch += 1 + if Lookup.Epoch == UInt32(0) + fill!(Lookup.Stamps, 0) + Lookup.Epoch = UInt32(1) + end end return nothing end @@ -60,7 +66,7 @@ end index = Int(HashCellIndex(Cell) & UInt(mask)) @inbounds for _ in 0:mask slot = index + 1 - if !Lookup.Filled[slot] + if Lookup.Stamps[slot] != Lookup.Epoch return default elseif Lookup.Keys[slot] == Cell return Lookup.Values[slot] @@ -75,8 +81,8 @@ end index = Int(HashCellIndex(Cell) & UInt(mask)) @inbounds for _ in 0:mask slot = index + 1 - if !Lookup.Filled[slot] - Lookup.Filled[slot] = true + if Lookup.Stamps[slot] != Lookup.Epoch + Lookup.Stamps[slot] = Lookup.Epoch Lookup.Keys[slot] = Cell Lookup.Values[slot] = value return nothing @@ -89,6 +95,26 @@ 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 @@ -163,16 +189,14 @@ function UpdateNeighbors!(Particles, InverseCutOff, ParticleRanges, Cells = @views Particles.Cells ParticleRanges[1] = 1 IndexCounter = 1 - ResetCellLookup!(CellLookup, length(Cells) * 2) + ResetCellLookup!(CellLookup, length(Cells) * 4) fill!(CellOffsets, zero(eltype(CellOffsets))) @inbounds for Index in eachindex(Cells) Cell = Cells[Index] - CellIndex = GetCellIndex(CellLookup, Cell, 0) - if CellIndex == 0 - IndexCounter += 1 - CellIndex = IndexCounter - SetCellIndex!(CellLookup, Cell, CellIndex) + CellIndex, IsNew = GetOrInsertCellIndex!(CellLookup, Cell, IndexCounter + 1) + if IsNew + IndexCounter = CellIndex UniqueCells[CellIndex] = Cell end CellOffsets[CellIndex] += 1 From deb8487458d4f267e8d4af0cb48a0d6a65ae51fc Mon Sep 17 00:00:00 2001 From: AhmedSalih3d <36305327+AhmedSalih3d@users.noreply.github.com> Date: Thu, 15 Jan 2026 21:23:40 +0100 Subject: [PATCH 4/6] Optimize neighbor list buffers --- src/SPHCellList.jl | 99 ++++++++++++++++++++++++++++-------------- src/SPHNeighborList.jl | 44 +++++++++++-------- 2 files changed, 92 insertions(+), 51 deletions(-) diff --git a/src/SPHCellList.jl b/src/SPHCellList.jl index 89111a3d..f6538d14 100644 --- a/src/SPHCellList.jl +++ b/src/SPHCellList.jl @@ -35,7 +35,7 @@ using LinearAlgebra function NeighborLoopPerParticle!(SimDensityDiffusion::SDD, SimViscosity::SV, SimKernel, SimMetaData::SimulationMetaData{D,T,NoShifting,NoKernelOutput,B,L}, SimConstants, SimParticles, ParticleRanges, - CellLookup, NeighborCellLists, dρdtI, + CellLookup, NeighborCellOffsets, NeighborCellCounts, NeighborCells, dρdtI, Acceleration, ∇Cᵢ, ∇◌rᵢ, max_visc = nothing, min_dt_force = nothing; @@ -55,8 +55,6 @@ using LinearAlgebra CellListIndex = GetCellIndex(CellLookup, CellIndex, 1) 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 +65,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 +92,7 @@ using LinearAlgebra function NeighborLoopPerParticle!(SimDensityDiffusion::SDD, SimViscosity::SV, SimKernel, SimMetaData::SimulationMetaData{D,T,NoShifting,K,B,L}, SimConstants, SimParticles, ParticleRanges, - CellLookup, NeighborCellLists, dρdtI, + CellLookup, NeighborCellOffsets, NeighborCellCounts, NeighborCells, dρdtI, Acceleration, ∇Cᵢ, ∇◌rᵢ, max_visc = nothing, min_dt_force = nothing; @@ -114,8 +115,6 @@ using LinearAlgebra CellListIndex = GetCellIndex(CellLookup, CellIndex, 1) 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 +127,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 +158,7 @@ using LinearAlgebra function NeighborLoopPerParticle!(SimDensityDiffusion::SDD, SimViscosity::SV, SimKernel, SimMetaData::SimulationMetaData{D,T,S,NoKernelOutput,B,L}, SimConstants, SimParticles, ParticleRanges, - CellLookup, NeighborCellLists, dρdtI, + CellLookup, NeighborCellOffsets, NeighborCellCounts, NeighborCells, dρdtI, Acceleration, ∇Cᵢ, ∇◌rᵢ, max_visc = nothing, min_dt_force = nothing; @@ -178,8 +180,6 @@ using LinearAlgebra CellListIndex = GetCellIndex(CellLookup, CellIndex, 1) 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, - CellLookup, NeighborCellLists, dρdtI, + CellLookup, NeighborCellOffsets, NeighborCellCounts, NeighborCells, dρdtI, Acceleration, ∇Cᵢ, ∇◌rᵢ, max_visc = nothing, min_dt_force = nothing; @@ -246,8 +249,6 @@ using LinearAlgebra CellListIndex = GetCellIndex(CellLookup, CellIndex, 1) 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 +261,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_ @@ -698,9 +702,9 @@ using LinearAlgebra SimMetaData::SimulationMetaData{Dimensions, FloatType, SMode, KMode, BMode, LMode}, SimConstants, SimParticles, FullStencil, ParticleRanges, UniqueCells, CellLookup, - ParticleOrder, CellOffsets, - NeighborCellLists, dρdtI, Velocityₙ⁺, - Positionₙ⁺, ρₙ⁺, ∇Cᵢ, ∇◌rᵢ, + ParticleOrder, CellOffsets, CellIndices, + NeighborCellOffsets, NeighborCellCounts, NeighborCells, + dρdtI, Velocityₙ⁺, Positionₙ⁺, ρₙ⁺, ∇Cᵢ, ∇◌rᵢ, MotionDefinition::Union{ Nothing, AbstractVector{ @@ -747,10 +751,10 @@ 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, CellLookup, ParticleOrder, CellOffsets) + @timeit SimMetaData.HourGlass "01a Actual Calculate IndexCounter" SimMetaData.IndexCounter = UpdateNeighbors!(SimParticles, SimKernel.H⁻¹, ParticleRanges, UniqueCells, CellLookup, ParticleOrder, CellOffsets, CellIndices) SimMetaData.Δx = zero(eltype(dρdtI)) UniqueCellsView = view(UniqueCells, 1:SimMetaData.IndexCounter) - BuildNeighborCellLists!(NeighborCellLists, FullStencil, UniqueCellsView, ParticleRanges, CellLookup) + BuildNeighborCellLists!(NeighborCellOffsets, NeighborCellCounts, NeighborCells, FullStencil, UniqueCellsView, ParticleRanges, CellLookup) end end @@ -767,14 +771,16 @@ using LinearAlgebra @timeit SimMetaData.HourGlass "04 First NeighborLoop" NeighborLoopPerParticle!( SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData, SimConstants, SimParticles, ParticleRanges, CellLookup, - NeighborCellLists, dρdtI, Acceleration, ∇Cᵢ, ∇◌rᵢ, + 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, CellLookup, - NeighborCellLists, dρdtI, Acceleration, ∇Cᵢ, ∇◌rᵢ, + NeighborCellOffsets, NeighborCellCounts, NeighborCells, dρdtI, + Acceleration, ∇Cᵢ, ∇◌rᵢ, ParticleOrder = ParticleOrder, ) end @@ -792,7 +798,8 @@ using LinearAlgebra @timeit SimMetaData.HourGlass "08 Second NeighborLoop" NeighborLoopPerParticle!( SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData, SimConstants, SimParticles, ParticleRanges, CellLookup, - NeighborCellLists, dρdtI, Acceleration, ∇Cᵢ, ∇◌rᵢ, + NeighborCellOffsets, NeighborCellCounts, NeighborCells, dρdtI, + Acceleration, ∇Cᵢ, ∇◌rᵢ, max_visc, min_dt_force, ParticleOrder = ParticleOrder, Position = Positionₙ⁺, @@ -803,7 +810,8 @@ using LinearAlgebra @timeit SimMetaData.HourGlass "08 Second NeighborLoop" NeighborLoopPerParticle!( SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData, SimConstants, SimParticles, ParticleRanges, CellLookup, - NeighborCellLists, dρdtI, Acceleration, ∇Cᵢ, ∇◌rᵢ, + NeighborCellOffsets, NeighborCellCounts, NeighborCells, dρdtI, + Acceleration, ∇Cᵢ, ∇◌rᵢ, max_visc, min_dt_force, ParticleOrder = ParticleOrder, Position = Positionₙ⁺, @@ -828,9 +836,17 @@ using LinearAlgebra end # To control if grid is exported in output - function DetermineOutput(::Val{true}, SimMetaData, output, ParticleRanges, NeighborCellLists, UniqueCellsView) + function DetermineOutput(::Val{true}, SimMetaData, output, ParticleRanges, + NeighborCellOffsets, NeighborCellCounts, NeighborCells, + UniqueCellsView) cell_particle_counts = ComputeCellParticleCounts(ParticleRanges, length(UniqueCellsView)) - cell_neighbor_counts = ComputeCellNeighborCounts(ParticleRanges, NeighborCellLists, length(UniqueCellsView)) + cell_neighbor_counts = ComputeCellNeighborCounts( + ParticleRanges, + NeighborCellOffsets, + NeighborCellCounts, + NeighborCells, + length(UniqueCellsView), + ) @timeit SimMetaData.HourGlass "13 Save Particle Data" begin output.enqueue_particles(SimMetaData.OutputIterationCounter) @@ -838,7 +854,9 @@ using LinearAlgebra end end - function DetermineOutput(::Val{false}, SimMetaData, output, ParticleRanges, NeighborCellLists, UniqueCellsView) + function DetermineOutput(::Val{false}, SimMetaData, output, ParticleRanges, + NeighborCellOffsets, NeighborCellCounts, NeighborCells, + UniqueCellsView) output.enqueue_particles(SimMetaData.OutputIterationCounter) end @@ -869,9 +887,12 @@ using LinearAlgebra UniqueCells = zeros(CartesianIndex{Dimensions}, NumberOfPoints) CellLookup = InitializeCellLookup(Val(Dimensions), NumberOfPoints * 2) 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)) ParticleOrder = zeros(Int, NumberOfPoints) CellOffsets = zeros(Int, length(ParticleRanges)) + CellIndices = zeros(Int, NumberOfPoints) output = SetupVTKOutput(SimMetaData, SimParticles, SimKernel, Dimensions) @@ -889,7 +910,9 @@ using LinearAlgebra ) cell_neighbor_counts = ComputeCellNeighborCounts( ParticleRanges, - NeighborCellLists, + NeighborCellOffsets, + NeighborCellCounts, + NeighborCells, SimMetaData.IndexCounter, ) end @@ -909,8 +932,9 @@ using LinearAlgebra @timeit SimMetaData.HourGlass "00 SimulationLoop" SimulationLoop( SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData, SimConstants, SimParticles, FullStencil, ParticleRanges, - UniqueCells, CellLookup, ParticleOrder, CellOffsets, - NeighborCellLists, dρdtI, Velocityₙ⁺, Positionₙ⁺, ρₙ⁺, + UniqueCells, CellLookup, ParticleOrder, CellOffsets, CellIndices, + NeighborCellOffsets, NeighborCellCounts, NeighborCells, + dρdtI, Velocityₙ⁺, Positionₙ⁺, ρₙ⁺, ∇Cᵢ, ∇◌rᵢ, MotionDefinition, ) push!(SimMetaData.TimeSteps, SimMetaData.CurrentTimeStep) @@ -921,7 +945,16 @@ 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, + UniqueCellsView, + ) if SimMetaData.TotalTime > SimMetaData.SimulationTime diff --git a/src/SPHNeighborList.jl b/src/SPHNeighborList.jl index 8fbb9920..cb22d325 100644 --- a/src/SPHNeighborList.jl +++ b/src/SPHNeighborList.jl @@ -119,22 +119,22 @@ function ConstructStencil(V::Val{d}) where d return CartesianIndices(ntuple(_ -> -1:1, V)) end -function BuildNeighborCellLists!(NeighborCellLists, FullStencil, UniqueCellsView, ParticleRanges, CellLookup) - TargetLen = length(UniqueCellsView) - OriginalLen = length(NeighborCellLists) - resize!(NeighborCellLists, TargetLen) - - if TargetLen > OriginalLen - @inbounds for Index in (OriginalLen + 1):TargetLen - NeighborCellLists[Index] = Int[] - end +function BuildNeighborCellLists!(NeighborCellOffsets, NeighborCellCounts, NeighborCells, + FullStencil, UniqueCellsView, ParticleRanges, CellLookup) + TargetLen = length(UniqueCellsView) + resize!(NeighborCellOffsets, TargetLen) + resize!(NeighborCellCounts, TargetLen) + + 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 = GetCellIndex(CellLookup, NeighborCell, 0) @@ -142,10 +142,13 @@ function BuildNeighborCellLists!(NeighborCellLists, FullStencil, UniqueCellsView StartIndex = ParticleRanges[NeighborIndex] EndIndex = ParticleRanges[NeighborIndex + 1] - 1 if StartIndex <= EndIndex - push!(Neighbors, NeighborIndex) + count += 1 + NeighborCells[write_index + count - 1] = NeighborIndex end end end + NeighborCellOffsets[CellIndex] = write_index + NeighborCellCounts[CellIndex] = count end return nothing @@ -183,13 +186,14 @@ 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, CellLookup, ParticleOrder, CellOffsets) + UniqueCells, CellLookup, ParticleOrder, CellOffsets, + CellIndices) ExtractCells!(Particles, InverseCutOff) Cells = @views Particles.Cells ParticleRanges[1] = 1 IndexCounter = 1 - ResetCellLookup!(CellLookup, length(Cells) * 4) + ResetCellLookup!(CellLookup, length(Cells) * 8) fill!(CellOffsets, zero(eltype(CellOffsets))) @inbounds for Index in eachindex(Cells) @@ -199,6 +203,7 @@ function UpdateNeighbors!(Particles, InverseCutOff, ParticleRanges, IndexCounter = CellIndex UniqueCells[CellIndex] = Cell end + CellIndices[Index] = CellIndex CellOffsets[CellIndex] += 1 end @@ -212,7 +217,7 @@ function UpdateNeighbors!(Particles, InverseCutOff, ParticleRanges, ParticleRanges[IndexCounter + 1] = RunningIndex @inbounds for Index in eachindex(Cells) - CellIndex = GetCellIndex(CellLookup, Cells[Index], 0) + CellIndex = CellIndices[Index] TargetIndex = CellOffsets[CellIndex] ParticleOrder[TargetIndex] = Index CellOffsets[CellIndex] = TargetIndex + 1 @@ -229,12 +234,15 @@ function ComputeCellParticleCounts(ParticleRanges, CellCount) return Counts end -function ComputeCellNeighborCounts(ParticleRanges, NeighborCellLists, CellCount) +function ComputeCellNeighborCounts(ParticleRanges, NeighborCellOffsets, NeighborCellCounts, NeighborCells, CellCount) Counts = ComputeCellParticleCounts(ParticleRanges, CellCount) Neighbors = Vector{Int}(undef, 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 From f2043d2cfd1b8199e21d5aeac5f99a253fbf50bd Mon Sep 17 00:00:00 2001 From: AhmedSalih3d <36305327+AhmedSalih3d@users.noreply.github.com> Date: Thu, 15 Jan 2026 21:49:38 +0100 Subject: [PATCH 5/6] Add dense grid neighbor lookup --- src/SPHCellList.jl | 74 +++++++++++++++++++++++++-------------- src/SPHNeighborList.jl | 79 +++++++++++++++++++++++++++++++++++++++--- 2 files changed, 122 insertions(+), 31 deletions(-) diff --git a/src/SPHCellList.jl b/src/SPHCellList.jl index f6538d14..7804c612 100644 --- a/src/SPHCellList.jl +++ b/src/SPHCellList.jl @@ -18,7 +18,9 @@ using ..OpenExternalPrograms using ..SPHKernels using ..SPHViscosityModels using ..SPHDensityDiffusionModels -using ..SPHNeighborList: BuildNeighborCellLists!, CellLookup, ComputeCellNeighborCounts, ComputeCellParticleCounts, ConstructStencil, ExtractCells!, GetCellIndex, InitializeCellLookup, MapFloor, UpdateNeighbors!, UpdateΔx! +using ..SPHNeighborList: BuildNeighborCellLists!, CellLookup, DenseCellLookup, ComputeCellNeighborCounts, + ComputeCellParticleCounts, ConstructStencil, ExtractCells!, GetCellIndexDispatch, + InitializeCellLookup, InitializeDenseCellLookup, MapFloor, UpdateNeighbors!, UpdateΔx! using StaticArrays using StructArrays: StructArray, foreachfield @@ -35,7 +37,7 @@ using LinearAlgebra function NeighborLoopPerParticle!(SimDensityDiffusion::SDD, SimViscosity::SV, SimKernel, SimMetaData::SimulationMetaData{D,T,NoShifting,NoKernelOutput,B,L}, SimConstants, SimParticles, ParticleRanges, - CellLookup, NeighborCellOffsets, NeighborCellCounts, NeighborCells, dρdtI, + CellIndices, NeighborCellOffsets, NeighborCellCounts, NeighborCells, dρdtI, Acceleration, ∇Cᵢ, ∇◌rᵢ, max_visc = nothing, min_dt_force = nothing; @@ -51,8 +53,7 @@ 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 = GetCellIndex(CellLookup, CellIndex, 1) + CellListIndex = CellIndices[i] SameCellStart = ParticleRanges[CellListIndex] SameCellEnd = ParticleRanges[CellListIndex + 1] - 1 @inbounds for j in SameCellStart:SameCellEnd @@ -92,7 +93,7 @@ using LinearAlgebra function NeighborLoopPerParticle!(SimDensityDiffusion::SDD, SimViscosity::SV, SimKernel, SimMetaData::SimulationMetaData{D,T,NoShifting,K,B,L}, SimConstants, SimParticles, ParticleRanges, - CellLookup, NeighborCellOffsets, NeighborCellCounts, NeighborCells, dρdtI, + CellIndices, NeighborCellOffsets, NeighborCellCounts, NeighborCells, dρdtI, Acceleration, ∇Cᵢ, ∇◌rᵢ, max_visc = nothing, min_dt_force = nothing; @@ -111,8 +112,7 @@ using LinearAlgebra acc_acc = zero(Acceleration[i]) kernel_acc = zero(Kernel[i]) kernel_grad_acc = zero(KernelGradient[i]) - CellIndex = Cells[i] - CellListIndex = GetCellIndex(CellLookup, CellIndex, 1) + CellListIndex = CellIndices[i] SameCellStart = ParticleRanges[CellListIndex] SameCellEnd = ParticleRanges[CellListIndex + 1] - 1 @inbounds for j in SameCellStart:SameCellEnd @@ -158,7 +158,7 @@ using LinearAlgebra function NeighborLoopPerParticle!(SimDensityDiffusion::SDD, SimViscosity::SV, SimKernel, SimMetaData::SimulationMetaData{D,T,S,NoKernelOutput,B,L}, SimConstants, SimParticles, ParticleRanges, - CellLookup, NeighborCellOffsets, NeighborCellCounts, NeighborCells, dρdtI, + CellIndices, NeighborCellOffsets, NeighborCellCounts, NeighborCells, dρdtI, Acceleration, ∇Cᵢ, ∇◌rᵢ, max_visc = nothing, min_dt_force = nothing; @@ -176,8 +176,7 @@ using LinearAlgebra acc_acc = zero(Acceleration[i]) shift_c_acc = zero(∇Cᵢ[i]) shift_r_acc = zero(∇◌rᵢ[i]) - CellIndex = Cells[i] - CellListIndex = GetCellIndex(CellLookup, CellIndex, 1) + CellListIndex = CellIndices[i] SameCellStart = ParticleRanges[CellListIndex] SameCellEnd = ParticleRanges[CellListIndex + 1] - 1 @inbounds for j in SameCellStart:SameCellEnd @@ -223,7 +222,7 @@ using LinearAlgebra function NeighborLoopPerParticle!(SimDensityDiffusion::SDD, SimViscosity::SV, SimKernel, SimMetaData::SimulationMetaData{D,T,S,K,B,L}, SimConstants, SimParticles, ParticleRanges, - CellLookup, NeighborCellOffsets, NeighborCellCounts, NeighborCells, dρdtI, + CellIndices, NeighborCellOffsets, NeighborCellCounts, NeighborCells, dρdtI, Acceleration, ∇Cᵢ, ∇◌rᵢ, max_visc = nothing, min_dt_force = nothing; @@ -245,8 +244,7 @@ using LinearAlgebra kernel_grad_acc = zero(KernelGradient[i]) shift_c_acc = zero(∇Cᵢ[i]) shift_r_acc = zero(∇◌rᵢ[i]) - CellIndex = Cells[i] - CellListIndex = GetCellIndex(CellLookup, CellIndex, 1) + CellListIndex = CellIndices[i] SameCellStart = ParticleRanges[CellListIndex] SameCellEnd = ParticleRanges[CellListIndex + 1] - 1 @inbounds for j in SameCellStart:SameCellEnd @@ -295,7 +293,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, CellLookup, Position, + SimConstants, ParticleRanges, CellLookup, DenseLookup, UseDense, Position, Density, GhostPoints, GhostNormals, ParticleType, bᵧ, Aᵧ; ParticleOrder) where {Dimensions, FloatType, SMode, KMode, BMode, LMode} @@ -318,7 +316,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 = GetCellIndex(CellLookup, SCellIndex, 1) + NeighborIdx = GetCellIndexDispatch(UseDense, DenseLookup, CellLookup, SCellIndex, 1) StartIndex_ = ParticleRanges[NeighborIdx] EndIndex_ = ParticleRanges[NeighborIdx + 1] - 1 @@ -636,7 +634,7 @@ using LinearAlgebra function ApplyMDBCBeforeHalf!(SimMetaData::SimulationMetaData{D,T,S,K,SimpleMDBC,L}, SimKernel, SimConstants, SimParticles, - ParticleRanges, CellLookup, Position, Density, + ParticleRanges, CellLookup, DenseLookup, UseDense, Position, Density, GhostPoints, GhostNormals, ParticleType; ParticleOrder ) where {D,T,S<:ShiftingMode,K<:KernelOutputMode,L<:LogMode} @@ -644,7 +642,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, CellLookup, 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 @@ -701,7 +699,7 @@ using LinearAlgebra @inbounds function SimulationLoop(SimDensityDiffusion::SDD, SimViscosity::SV, SimKernel, SimMetaData::SimulationMetaData{Dimensions, FloatType, SMode, KMode, BMode, LMode}, SimConstants, SimParticles, FullStencil, - ParticleRanges, UniqueCells, CellLookup, + ParticleRanges, UniqueCells, CellLookup, DenseLookup, UseDense, ParticleOrder, CellOffsets, CellIndices, NeighborCellOffsets, NeighborCellCounts, NeighborCells, dρdtI, Velocityₙ⁺, Positionₙ⁺, ρₙ⁺, ∇Cᵢ, ∇◌rᵢ, @@ -751,10 +749,32 @@ 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, CellLookup, ParticleOrder, CellOffsets, CellIndices) + @timeit SimMetaData.HourGlass "01a Actual Calculate IndexCounter" begin + SimMetaData.IndexCounter, UseDense[] = UpdateNeighbors!( + SimParticles, + SimKernel.H⁻¹, + ParticleRanges, + UniqueCells, + CellLookup, + ParticleOrder, + CellOffsets, + CellIndices, + DenseLookup, + ) + end SimMetaData.Δx = zero(eltype(dρdtI)) UniqueCellsView = view(UniqueCells, 1:SimMetaData.IndexCounter) - BuildNeighborCellLists!(NeighborCellOffsets, NeighborCellCounts, NeighborCells, FullStencil, UniqueCellsView, ParticleRanges, CellLookup) + BuildNeighborCellLists!( + NeighborCellOffsets, + NeighborCellCounts, + NeighborCells, + FullStencil, + UniqueCellsView, + ParticleRanges, + CellLookup, + DenseLookup, + Val(UseDense[]), + ) end end @@ -764,13 +784,13 @@ 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, CellLookup, 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, CellLookup, + SimConstants, SimParticles, ParticleRanges, CellIndices, NeighborCellOffsets, NeighborCellCounts, NeighborCells, dρdtI, Acceleration, ∇Cᵢ, ∇◌rᵢ, ParticleOrder = ParticleOrder, @@ -778,7 +798,7 @@ using LinearAlgebra else @timeit SimMetaData.HourGlass "04 First NeighborLoop" NeighborLoopPerParticle!( SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData, - SimConstants, SimParticles, ParticleRanges, CellLookup, + SimConstants, SimParticles, ParticleRanges, CellIndices, NeighborCellOffsets, NeighborCellCounts, NeighborCells, dρdtI, Acceleration, ∇Cᵢ, ∇◌rᵢ, ParticleOrder = ParticleOrder, @@ -797,7 +817,7 @@ 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, CellLookup, + SimConstants, SimParticles, ParticleRanges, CellIndices, NeighborCellOffsets, NeighborCellCounts, NeighborCells, dρdtI, Acceleration, ∇Cᵢ, ∇◌rᵢ, max_visc, min_dt_force, @@ -809,7 +829,7 @@ using LinearAlgebra else @timeit SimMetaData.HourGlass "08 Second NeighborLoop" NeighborLoopPerParticle!( SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData, - SimConstants, SimParticles, ParticleRanges, CellLookup, + SimConstants, SimParticles, ParticleRanges, CellIndices, NeighborCellOffsets, NeighborCellCounts, NeighborCells, dρdtI, Acceleration, ∇Cᵢ, ∇◌rᵢ, max_visc, min_dt_force, @@ -886,6 +906,8 @@ using LinearAlgebra ParticleRanges = zeros(Int, NumberOfPoints + 1 + 1) # +1 for the last particle, +1 for dummy entry UniqueCells = zeros(CartesianIndex{Dimensions}, NumberOfPoints) CellLookup = InitializeCellLookup(Val(Dimensions), NumberOfPoints * 2) + DenseLookup = InitializeDenseCellLookup(Val(Dimensions)) + UseDense = Ref(false) FullStencil = ConstructStencil(Val(Dimensions)) NeighborCellOffsets = zeros(Int, length(UniqueCells)) NeighborCellCounts = zeros(Int, length(UniqueCells)) @@ -932,7 +954,7 @@ using LinearAlgebra @timeit SimMetaData.HourGlass "00 SimulationLoop" SimulationLoop( SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData, SimConstants, SimParticles, FullStencil, ParticleRanges, - UniqueCells, CellLookup, ParticleOrder, CellOffsets, CellIndices, + UniqueCells, CellLookup, DenseLookup, UseDense, ParticleOrder, CellOffsets, CellIndices, NeighborCellOffsets, NeighborCellCounts, NeighborCells, dρdtI, Velocityₙ⁺, Positionₙ⁺, ρₙ⁺, ∇Cᵢ, ∇◌rᵢ, MotionDefinition, diff --git a/src/SPHNeighborList.jl b/src/SPHNeighborList.jl index cb22d325..1e9ee57b 100644 --- a/src/SPHNeighborList.jl +++ b/src/SPHNeighborList.jl @@ -2,7 +2,8 @@ module SPHNeighborList export ConstructStencil, ExtractCells!, UpdateNeighbors!, BuildNeighborCellLists!, ComputeCellParticleCounts, ComputeCellNeighborCounts, UpdateΔx!, - CellLookup, InitializeCellLookup, ResetCellLookup!, GetCellIndex, SetCellIndex! + CellLookup, DenseCellLookup, InitializeCellLookup, InitializeDenseCellLookup, + ResetCellLookup!, GetCellIndex, SetCellIndex!, GetCellIndexDense, GetCellIndexDispatch using StaticArrays @@ -14,6 +15,13 @@ mutable struct CellLookup{D} Epoch::UInt32 end +mutable struct DenseCellLookup{D} + MinCell::SVector{D, Int} + Dims::SVector{D, Int} + Strides::SVector{D, Int} + IndexMap::Vector{Int} +end + @inline function ZeroCell(::Val{D}) where D return CartesianIndex(ntuple(_ -> 0, Val(D))) end @@ -32,6 +40,41 @@ function InitializeCellLookup(::Val{D}, capacity::Int) where 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 @@ -120,7 +163,8 @@ function ConstructStencil(V::Val{d}) where d end function BuildNeighborCellLists!(NeighborCellOffsets, NeighborCellCounts, NeighborCells, - FullStencil, UniqueCellsView, ParticleRanges, CellLookup) + FullStencil, UniqueCellsView, ParticleRanges, + CellLookup, DenseLookup, UseDense) TargetLen = length(UniqueCellsView) resize!(NeighborCellOffsets, TargetLen) resize!(NeighborCellCounts, TargetLen) @@ -137,7 +181,7 @@ function BuildNeighborCellLists!(NeighborCellOffsets, NeighborCellCounts, Neighb count = 0 for Offset in FullStencil NeighborCell = Cell + Offset - NeighborIndex = GetCellIndex(CellLookup, NeighborCell, 0) + NeighborIndex = GetCellIndexDispatch(UseDense, DenseLookup, CellLookup, NeighborCell, 0) if NeighborIndex != 0 && NeighborIndex != CellIndex StartIndex = ParticleRanges[NeighborIndex] EndIndex = ParticleRanges[NeighborIndex + 1] - 1 @@ -187,7 +231,7 @@ without reordering the particle arrays. """ function UpdateNeighbors!(Particles, InverseCutOff, ParticleRanges, UniqueCells, CellLookup, ParticleOrder, CellOffsets, - CellIndices) + CellIndices, DenseLookup) ExtractCells!(Particles, InverseCutOff) Cells = @views Particles.Cells @@ -195,9 +239,15 @@ function UpdateNeighbors!(Particles, InverseCutOff, ParticleRanges, IndexCounter = 1 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] + 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 @@ -223,7 +273,26 @@ function UpdateNeighbors!(Particles, InverseCutOff, ParticleRanges, 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) From e2d19cf5a27e411c0a5822bf837fbd83e68b4128 Mon Sep 17 00:00:00 2001 From: AhmedSalih3d <36305327+AhmedSalih3d@users.noreply.github.com> Date: Thu, 15 Jan 2026 22:04:19 +0100 Subject: [PATCH 6/6] Add neighbor list scratch buffers --- src/SPHCellList.jl | 57 ++++++++++++++++++++++++++---------------- src/SPHExample.jl | 4 ++- src/SPHNeighborList.jl | 56 ++++++++++++++++++++++++++++++++++------- 3 files changed, 85 insertions(+), 32 deletions(-) diff --git a/src/SPHCellList.jl b/src/SPHCellList.jl index 7804c612..760678e4 100644 --- a/src/SPHCellList.jl +++ b/src/SPHCellList.jl @@ -18,9 +18,10 @@ using ..OpenExternalPrograms using ..SPHKernels using ..SPHViscosityModels using ..SPHDensityDiffusionModels -using ..SPHNeighborList: BuildNeighborCellLists!, CellLookup, DenseCellLookup, ComputeCellNeighborCounts, - ComputeCellParticleCounts, ConstructStencil, ExtractCells!, GetCellIndexDispatch, - InitializeCellLookup, InitializeDenseCellLookup, 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 @@ -702,7 +703,7 @@ using LinearAlgebra ParticleRanges, UniqueCells, CellLookup, DenseLookup, UseDense, ParticleOrder, CellOffsets, CellIndices, NeighborCellOffsets, NeighborCellCounts, NeighborCells, - dρdtI, Velocityₙ⁺, Positionₙ⁺, ρₙ⁺, ∇Cᵢ, ∇◌rᵢ, + NeighborScratch, dρdtI, Velocityₙ⁺, Positionₙ⁺, ρₙ⁺, ∇Cᵢ, ∇◌rᵢ, MotionDefinition::Union{ Nothing, AbstractVector{ @@ -750,18 +751,19 @@ using LinearAlgebra # 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" begin - SimMetaData.IndexCounter, UseDense[] = UpdateNeighbors!( - SimParticles, - SimKernel.H⁻¹, - ParticleRanges, + SimMetaData.IndexCounter, UseDense[] = UpdateNeighbors!( + SimParticles, + SimKernel.H⁻¹, + ParticleRanges, UniqueCells, CellLookup, ParticleOrder, CellOffsets, - CellIndices, - DenseLookup, - ) + CellIndices, + DenseLookup, + ) end + UpdateCellCounts!(NeighborScratch, ParticleRanges, SimMetaData.IndexCounter) SimMetaData.Δx = zero(eltype(dρdtI)) UniqueCellsView = view(UniqueCells, 1:SimMetaData.IndexCounter) BuildNeighborCellLists!( @@ -774,6 +776,7 @@ using LinearAlgebra CellLookup, DenseLookup, Val(UseDense[]), + NeighborScratch.CellOccupied, ) end end @@ -858,10 +861,12 @@ using LinearAlgebra # To control if grid is exported in output function DetermineOutput(::Val{true}, SimMetaData, output, ParticleRanges, NeighborCellOffsets, NeighborCellCounts, NeighborCells, + NeighborScratch, UniqueCellsView) - cell_particle_counts = ComputeCellParticleCounts(ParticleRanges, length(UniqueCellsView)) - cell_neighbor_counts = ComputeCellNeighborCounts( - ParticleRanges, + UpdateCellCounts!(NeighborScratch, ParticleRanges, length(UniqueCellsView)) + ComputeCellNeighborCounts!( + NeighborScratch.NeighborCounts, + NeighborScratch.CellCounts, NeighborCellOffsets, NeighborCellCounts, NeighborCells, @@ -870,12 +875,18 @@ using LinearAlgebra @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, NeighborCellOffsets, NeighborCellCounts, NeighborCells, + NeighborScratch, UniqueCellsView) output.enqueue_particles(SimMetaData.OutputIterationCounter) end @@ -912,6 +923,7 @@ using LinearAlgebra 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) @@ -926,17 +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, + 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, @@ -955,7 +967,7 @@ using LinearAlgebra SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData, SimConstants, SimParticles, FullStencil, ParticleRanges, UniqueCells, CellLookup, DenseLookup, UseDense, ParticleOrder, CellOffsets, CellIndices, - NeighborCellOffsets, NeighborCellCounts, NeighborCells, + NeighborCellOffsets, NeighborCellCounts, NeighborCells, NeighborScratch, dρdtI, Velocityₙ⁺, Positionₙ⁺, ρₙ⁺, ∇Cᵢ, ∇◌rᵢ, MotionDefinition, ) @@ -975,6 +987,7 @@ using LinearAlgebra NeighborCellOffsets, NeighborCellCounts, NeighborCells, + NeighborScratch, UniqueCellsView, ) 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 1e9ee57b..64bce3ff 100644 --- a/src/SPHNeighborList.jl +++ b/src/SPHNeighborList.jl @@ -3,7 +3,8 @@ module SPHNeighborList export ConstructStencil, ExtractCells!, UpdateNeighbors!, BuildNeighborCellLists!, ComputeCellParticleCounts, ComputeCellNeighborCounts, UpdateΔx!, CellLookup, DenseCellLookup, InitializeCellLookup, InitializeDenseCellLookup, - ResetCellLookup!, GetCellIndex, SetCellIndex!, GetCellIndexDense, GetCellIndexDispatch + ResetCellLookup!, GetCellIndex, SetCellIndex!, GetCellIndexDense, GetCellIndexDispatch, + NeighborListScratch, UpdateCellCounts!, ComputeCellParticleCounts!, ComputeCellNeighborCounts! using StaticArrays @@ -22,6 +23,23 @@ mutable struct DenseCellLookup{D} 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 @@ -164,7 +182,7 @@ end function BuildNeighborCellLists!(NeighborCellOffsets, NeighborCellCounts, NeighborCells, FullStencil, UniqueCellsView, ParticleRanges, - CellLookup, DenseLookup, UseDense) + CellLookup, DenseLookup, UseDense, CellOccupied) TargetLen = length(UniqueCellsView) resize!(NeighborCellOffsets, TargetLen) resize!(NeighborCellCounts, TargetLen) @@ -182,13 +200,9 @@ function BuildNeighborCellLists!(NeighborCellOffsets, NeighborCellCounts, Neighb for Offset in FullStencil NeighborCell = Cell + Offset NeighborIndex = GetCellIndexDispatch(UseDense, DenseLookup, CellLookup, NeighborCell, 0) - if NeighborIndex != 0 && NeighborIndex != CellIndex - StartIndex = ParticleRanges[NeighborIndex] - EndIndex = ParticleRanges[NeighborIndex + 1] - 1 - if StartIndex <= EndIndex - count += 1 - NeighborCells[write_index + count - 1] = NeighborIndex - end + if NeighborIndex != 0 && NeighborIndex != CellIndex && CellOccupied[NeighborIndex] + count += 1 + NeighborCells[write_index + count - 1] = NeighborIndex end end NeighborCellOffsets[CellIndex] = write_index @@ -297,6 +311,12 @@ 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 @@ -306,6 +326,12 @@ end 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 count = NeighborCellCounts[Index] @@ -319,6 +345,18 @@ function ComputeCellNeighborCounts(ParticleRanges, NeighborCellOffsets, Neighbor 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)