From e582504a16e15cc874a28bb9f0a83d0f4041e967 Mon Sep 17 00:00:00 2001 From: AhmedSalih3d <36305327+AhmedSalih3d@users.noreply.github.com> Date: Wed, 14 Jan 2026 22:23:02 +0100 Subject: [PATCH 1/8] Add non-sorting neighbor list option --- README.md | 3 + src/SPHCellList.jl | 331 ++++++++++++++++++------- src/SPHExample.jl | 3 +- src/SPHNeighborList.jl | 50 +++- src/SimulationMetaDataConfiguration.jl | 1 + 5 files changed, 290 insertions(+), 98 deletions(-) diff --git a/README.md b/README.md index e1bf6bb4..fb3cc80c 100644 --- a/README.md +++ b/README.md @@ -88,6 +88,9 @@ To color exported cell grids by particle counts, set `ExportGridCellParticleCounts=true` in `SimulationMetaData`. This also adds a `ParticleNeighborsPerCell` array that includes each cell's particle count minus one plus the particles in its neighbor stencil. +If you want to avoid re-sorting particles each neighbor rebuild, set +`UseSortingNeighborList=false` to build a per-cell particle ordering buffer +without reordering particle storage. ## Help diff --git a/src/SPHCellList.jl b/src/SPHCellList.jl index fef4e636..feedc1b8 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!, ComputeCellNeighborCounts, ComputeCellParticleCounts, ConstructStencil, ExtractCells!, MapFloor, UpdateNeighbors!, UpdateNeighborsNoSort!, UpdateΔx! using StaticArrays using StructArrays: StructArray, foreachfield @@ -32,6 +32,14 @@ using Base.Threads using LinearAlgebra using Bumper + @inline function ResolveParticleIndex(::Nothing, Index) + return Index + end + + @inline function ResolveParticleIndex(ParticleOrder, Index) + return @inbounds ParticleOrder[Index] + end + function NeighborLoopPerParticle!(SimDensityDiffusion::SDD, SimViscosity::SV, SimKernel, SimMetaData::SimulationMetaData{D,T,NoShifting,NoKernelOutput,B,L}, SimConstants, SimParticles, ParticleRanges, @@ -42,7 +50,8 @@ using LinearAlgebra Position = SimParticles.Position, Density = SimParticles.Density, Pressure = SimParticles.Pressure, - Velocity = SimParticles.Velocity) where {D,T, + Velocity = SimParticles.Velocity, + ParticleOrder = nothing) where {D,T, B<:MDBCMode,L<:LogMode, SDD<:SPHDensityDiffusion, SV<:SPHViscosity} @@ -56,30 +65,55 @@ using LinearAlgebra SameCellEnd = ParticleRanges[CellListIndex + 1] - 1 NeighborCellIndices = NeighborCellLists[CellListIndex] - @inbounds for j in SameCellStart:(i - 1) - dρdt_acc, acc_acc = ComputeInteractionsPerParticleNoKernel!( - SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData, - SimConstants, SimParticles, Position, Density, Pressure, - Velocity, MotionLimiter, dρdt_acc, acc_acc, i, j, - ) - end - @inbounds for j in (i + 1):SameCellEnd - dρdt_acc, acc_acc = ComputeInteractionsPerParticleNoKernel!( - SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData, - SimConstants, SimParticles, Position, Density, Pressure, - Velocity, MotionLimiter, dρdt_acc, acc_acc, i, j, - ) - end - for NeighborIdx in NeighborCellIndices - StartIndex_ = ParticleRanges[NeighborIdx] - EndIndex_ = ParticleRanges[NeighborIdx + 1] - 1 - @inbounds for j in StartIndex_:EndIndex_ + if ParticleOrder === nothing + @inbounds for j in SameCellStart:(i - 1) + dρdt_acc, acc_acc = ComputeInteractionsPerParticleNoKernel!( + SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData, + SimConstants, SimParticles, Position, Density, Pressure, + Velocity, MotionLimiter, dρdt_acc, acc_acc, i, j, + ) + end + @inbounds for j in (i + 1):SameCellEnd dρdt_acc, acc_acc = ComputeInteractionsPerParticleNoKernel!( SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData, SimConstants, SimParticles, Position, Density, Pressure, Velocity, MotionLimiter, dρdt_acc, acc_acc, i, j, ) end + for NeighborIdx in NeighborCellIndices + StartIndex_ = ParticleRanges[NeighborIdx] + EndIndex_ = ParticleRanges[NeighborIdx + 1] - 1 + @inbounds for j in StartIndex_:EndIndex_ + dρdt_acc, acc_acc = ComputeInteractionsPerParticleNoKernel!( + SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData, + SimConstants, SimParticles, Position, Density, Pressure, + Velocity, MotionLimiter, dρdt_acc, acc_acc, i, j, + ) + end + end + else + @inbounds for j in SameCellStart:SameCellEnd + jIndex = ResolveParticleIndex(ParticleOrder, j) + if jIndex != i + dρdt_acc, acc_acc = ComputeInteractionsPerParticleNoKernel!( + SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData, + SimConstants, SimParticles, Position, Density, Pressure, + Velocity, MotionLimiter, dρdt_acc, acc_acc, i, jIndex, + ) + end + end + for NeighborIdx in NeighborCellIndices + StartIndex_ = ParticleRanges[NeighborIdx] + EndIndex_ = ParticleRanges[NeighborIdx + 1] - 1 + @inbounds for j in StartIndex_:EndIndex_ + jIndex = ResolveParticleIndex(ParticleOrder, j) + dρdt_acc, acc_acc = ComputeInteractionsPerParticleNoKernel!( + SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData, + SimConstants, SimParticles, Position, Density, Pressure, + Velocity, MotionLimiter, dρdt_acc, acc_acc, i, jIndex, + ) + end + end end dρdtI[i] = dρdt_acc @@ -100,7 +134,8 @@ using LinearAlgebra Position = SimParticles.Position, Density = SimParticles.Density, Pressure = SimParticles.Pressure, - Velocity = SimParticles.Velocity) where {D,T, + Velocity = SimParticles.Velocity, + ParticleOrder = nothing) where {D,T, K<:KernelOutputMode, B<:MDBCMode,L<:LogMode, SDD<:SPHDensityDiffusion, @@ -117,28 +152,8 @@ using LinearAlgebra SameCellEnd = ParticleRanges[CellListIndex + 1] - 1 NeighborCellIndices = NeighborCellLists[CellListIndex] - @inbounds for j in SameCellStart:(i - 1) - dρdt_acc, acc_acc, kernel_acc, kernel_grad_acc = - ComputeInteractionsPerParticle!( - SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData, - SimConstants, SimParticles, Position, Density, Pressure, - Velocity, MotionLimiter, dρdt_acc, acc_acc, kernel_acc, - kernel_grad_acc, i, j, - ) - end - @inbounds for j in (i + 1):SameCellEnd - dρdt_acc, acc_acc, kernel_acc, kernel_grad_acc = - ComputeInteractionsPerParticle!( - SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData, - SimConstants, SimParticles, Position, Density, Pressure, - Velocity, MotionLimiter, dρdt_acc, acc_acc, kernel_acc, - kernel_grad_acc, i, j, - ) - end - for NeighborIdx in NeighborCellIndices - StartIndex_ = ParticleRanges[NeighborIdx] - EndIndex_ = ParticleRanges[NeighborIdx + 1] - 1 - @inbounds for j in StartIndex_:EndIndex_ + if ParticleOrder === nothing + @inbounds for j in SameCellStart:(i - 1) dρdt_acc, acc_acc, kernel_acc, kernel_grad_acc = ComputeInteractionsPerParticle!( SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData, @@ -147,6 +162,55 @@ using LinearAlgebra kernel_grad_acc, i, j, ) end + @inbounds for j in (i + 1):SameCellEnd + dρdt_acc, acc_acc, kernel_acc, kernel_grad_acc = + ComputeInteractionsPerParticle!( + SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData, + SimConstants, SimParticles, Position, Density, Pressure, + Velocity, MotionLimiter, dρdt_acc, acc_acc, kernel_acc, + kernel_grad_acc, i, j, + ) + end + for NeighborIdx in NeighborCellIndices + StartIndex_ = ParticleRanges[NeighborIdx] + EndIndex_ = ParticleRanges[NeighborIdx + 1] - 1 + @inbounds for j in StartIndex_:EndIndex_ + dρdt_acc, acc_acc, kernel_acc, kernel_grad_acc = + ComputeInteractionsPerParticle!( + SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData, + SimConstants, SimParticles, Position, Density, Pressure, + Velocity, MotionLimiter, dρdt_acc, acc_acc, kernel_acc, + kernel_grad_acc, i, j, + ) + end + end + else + @inbounds for j in SameCellStart:SameCellEnd + jIndex = ResolveParticleIndex(ParticleOrder, j) + if jIndex != i + dρdt_acc, acc_acc, kernel_acc, kernel_grad_acc = + ComputeInteractionsPerParticle!( + SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData, + SimConstants, SimParticles, Position, Density, Pressure, + Velocity, MotionLimiter, dρdt_acc, acc_acc, kernel_acc, + kernel_grad_acc, i, jIndex, + ) + end + end + for NeighborIdx in NeighborCellIndices + StartIndex_ = ParticleRanges[NeighborIdx] + EndIndex_ = ParticleRanges[NeighborIdx + 1] - 1 + @inbounds for j in StartIndex_:EndIndex_ + jIndex = ResolveParticleIndex(ParticleOrder, j) + dρdt_acc, acc_acc, kernel_acc, kernel_grad_acc = + ComputeInteractionsPerParticle!( + SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData, + SimConstants, SimParticles, Position, Density, Pressure, + Velocity, MotionLimiter, dρdt_acc, acc_acc, kernel_acc, + kernel_grad_acc, i, jIndex, + ) + end + end end dρdtI[i] = dρdt_acc @@ -169,7 +233,8 @@ using LinearAlgebra Position = SimParticles.Position, Density = SimParticles.Density, Pressure = SimParticles.Pressure, - Velocity = SimParticles.Velocity) where {D,T, + Velocity = SimParticles.Velocity, + ParticleOrder = nothing) where {D,T, S<:ShiftingMode,B<:MDBCMode, L<:LogMode,SDD<:SPHDensityDiffusion, SV<:SPHViscosity} @@ -185,28 +250,8 @@ using LinearAlgebra SameCellEnd = ParticleRanges[CellListIndex + 1] - 1 NeighborCellIndices = NeighborCellLists[CellListIndex] - @inbounds for j in SameCellStart:(i - 1) - dρdt_acc, acc_acc, shift_c_acc, shift_r_acc = - ComputeInteractionsPerParticleNoKernel!( - SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData, - SimConstants, SimParticles, Position, Density, Pressure, - Velocity, MotionLimiter, dρdt_acc, acc_acc, shift_c_acc, - shift_r_acc, i, j, - ) - end - @inbounds for j in (i + 1):SameCellEnd - dρdt_acc, acc_acc, shift_c_acc, shift_r_acc = - ComputeInteractionsPerParticleNoKernel!( - SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData, - SimConstants, SimParticles, Position, Density, Pressure, - Velocity, MotionLimiter, dρdt_acc, acc_acc, shift_c_acc, - shift_r_acc, i, j, - ) - end - for NeighborIdx in NeighborCellIndices - StartIndex_ = ParticleRanges[NeighborIdx] - EndIndex_ = ParticleRanges[NeighborIdx + 1] - 1 - @inbounds for j in StartIndex_:EndIndex_ + if ParticleOrder === nothing + @inbounds for j in SameCellStart:(i - 1) dρdt_acc, acc_acc, shift_c_acc, shift_r_acc = ComputeInteractionsPerParticleNoKernel!( SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData, @@ -215,6 +260,55 @@ using LinearAlgebra shift_r_acc, i, j, ) end + @inbounds for j in (i + 1):SameCellEnd + dρdt_acc, acc_acc, shift_c_acc, shift_r_acc = + ComputeInteractionsPerParticleNoKernel!( + SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData, + SimConstants, SimParticles, Position, Density, Pressure, + Velocity, MotionLimiter, dρdt_acc, acc_acc, shift_c_acc, + shift_r_acc, i, j, + ) + end + for NeighborIdx in NeighborCellIndices + StartIndex_ = ParticleRanges[NeighborIdx] + EndIndex_ = ParticleRanges[NeighborIdx + 1] - 1 + @inbounds for j in StartIndex_:EndIndex_ + dρdt_acc, acc_acc, shift_c_acc, shift_r_acc = + ComputeInteractionsPerParticleNoKernel!( + SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData, + SimConstants, SimParticles, Position, Density, Pressure, + Velocity, MotionLimiter, dρdt_acc, acc_acc, shift_c_acc, + shift_r_acc, i, j, + ) + end + end + else + @inbounds for j in SameCellStart:SameCellEnd + jIndex = ResolveParticleIndex(ParticleOrder, j) + if jIndex != i + dρdt_acc, acc_acc, shift_c_acc, shift_r_acc = + ComputeInteractionsPerParticleNoKernel!( + SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData, + SimConstants, SimParticles, Position, Density, Pressure, + Velocity, MotionLimiter, dρdt_acc, acc_acc, shift_c_acc, + shift_r_acc, i, jIndex, + ) + end + end + for NeighborIdx in NeighborCellIndices + StartIndex_ = ParticleRanges[NeighborIdx] + EndIndex_ = ParticleRanges[NeighborIdx + 1] - 1 + @inbounds for j in StartIndex_:EndIndex_ + jIndex = ResolveParticleIndex(ParticleOrder, j) + dρdt_acc, acc_acc, shift_c_acc, shift_r_acc = + ComputeInteractionsPerParticleNoKernel!( + SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData, + SimConstants, SimParticles, Position, Density, Pressure, + Velocity, MotionLimiter, dρdt_acc, acc_acc, shift_c_acc, + shift_r_acc, i, jIndex, + ) + end + end end dρdtI[i] = dρdt_acc @@ -237,7 +331,8 @@ using LinearAlgebra Position = SimParticles.Position, Density = SimParticles.Density, Pressure = SimParticles.Pressure, - Velocity = SimParticles.Velocity) where {D,T, + Velocity = SimParticles.Velocity, + ParticleOrder = nothing) where {D,T, S<:ShiftingMode, K<:KernelOutputMode, B<:MDBCMode,L<:LogMode, @@ -257,28 +352,8 @@ using LinearAlgebra SameCellEnd = ParticleRanges[CellListIndex + 1] - 1 NeighborCellIndices = NeighborCellLists[CellListIndex] - @inbounds for j in SameCellStart:(i - 1) - dρdt_acc, acc_acc, kernel_acc, kernel_grad_acc, shift_c_acc, - shift_r_acc = ComputeInteractionsPerParticle!( - SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData, - SimConstants, SimParticles, Position, Density, Pressure, - Velocity, MotionLimiter, dρdt_acc, acc_acc, kernel_acc, - kernel_grad_acc, shift_c_acc, shift_r_acc, i, j, - ) - end - @inbounds for j in (i + 1):SameCellEnd - dρdt_acc, acc_acc, kernel_acc, kernel_grad_acc, shift_c_acc, - shift_r_acc = ComputeInteractionsPerParticle!( - SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData, - SimConstants, SimParticles, Position, Density, Pressure, - Velocity, MotionLimiter, dρdt_acc, acc_acc, kernel_acc, - kernel_grad_acc, shift_c_acc, shift_r_acc, i, j, - ) - end - for NeighborIdx in NeighborCellIndices - StartIndex_ = ParticleRanges[NeighborIdx] - EndIndex_ = ParticleRanges[NeighborIdx + 1] - 1 - @inbounds for j in StartIndex_:EndIndex_ + if ParticleOrder === nothing + @inbounds for j in SameCellStart:(i - 1) dρdt_acc, acc_acc, kernel_acc, kernel_grad_acc, shift_c_acc, shift_r_acc = ComputeInteractionsPerParticle!( SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData, @@ -287,6 +362,55 @@ using LinearAlgebra kernel_grad_acc, shift_c_acc, shift_r_acc, i, j, ) end + @inbounds for j in (i + 1):SameCellEnd + dρdt_acc, acc_acc, kernel_acc, kernel_grad_acc, shift_c_acc, + shift_r_acc = ComputeInteractionsPerParticle!( + SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData, + SimConstants, SimParticles, Position, Density, Pressure, + Velocity, MotionLimiter, dρdt_acc, acc_acc, kernel_acc, + kernel_grad_acc, shift_c_acc, shift_r_acc, i, j, + ) + end + for NeighborIdx in NeighborCellIndices + StartIndex_ = ParticleRanges[NeighborIdx] + EndIndex_ = ParticleRanges[NeighborIdx + 1] - 1 + @inbounds for j in StartIndex_:EndIndex_ + dρdt_acc, acc_acc, kernel_acc, kernel_grad_acc, shift_c_acc, + shift_r_acc = ComputeInteractionsPerParticle!( + SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData, + SimConstants, SimParticles, Position, Density, Pressure, + Velocity, MotionLimiter, dρdt_acc, acc_acc, kernel_acc, + kernel_grad_acc, shift_c_acc, shift_r_acc, i, j, + ) + end + end + else + @inbounds for j in SameCellStart:SameCellEnd + jIndex = ResolveParticleIndex(ParticleOrder, j) + if jIndex != i + dρdt_acc, acc_acc, kernel_acc, kernel_grad_acc, shift_c_acc, + shift_r_acc = ComputeInteractionsPerParticle!( + SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData, + SimConstants, SimParticles, Position, Density, Pressure, + Velocity, MotionLimiter, dρdt_acc, acc_acc, kernel_acc, + kernel_grad_acc, shift_c_acc, shift_r_acc, i, jIndex, + ) + end + end + for NeighborIdx in NeighborCellIndices + StartIndex_ = ParticleRanges[NeighborIdx] + EndIndex_ = ParticleRanges[NeighborIdx + 1] - 1 + @inbounds for j in StartIndex_:EndIndex_ + jIndex = ResolveParticleIndex(ParticleOrder, j) + dρdt_acc, acc_acc, kernel_acc, kernel_grad_acc, shift_c_acc, + shift_r_acc = ComputeInteractionsPerParticle!( + SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData, + SimConstants, SimParticles, Position, Density, Pressure, + Velocity, MotionLimiter, dρdt_acc, acc_acc, kernel_acc, + kernel_grad_acc, shift_c_acc, shift_r_acc, i, jIndex, + ) + end + end end dρdtI[i] = dρdt_acc @@ -307,7 +431,8 @@ using LinearAlgebra SimMetaData::SimulationMetaData{Dimensions, FloatType, SMode, KMode, BMode, LMode}, SimConstants, ParticleRanges, CellDict, Position, Density, GhostPoints, GhostNormals, ParticleType, - bᵧ, Aᵧ) where {Dimensions, FloatType, SMode, KMode, BMode, LMode} + bᵧ, Aᵧ; + ParticleOrder = nothing) where {Dimensions, FloatType, SMode, KMode, BMode, LMode} FullStencil = ConstructStencil(Val(Dimensions)) @@ -334,9 +459,10 @@ using LinearAlgebra for j in StartIndex_:EndIndex_ # change ComputeInteractions to take & return contributions, e.g.: + jIndex = ResolveParticleIndex(ParticleOrder, j) bΔ, AΔ = ComputeInteractionsMDBC!(SimKernel, SimMetaData, SimConstants, Position, Density, ParticleType, - GhostPoints, iter, j) + GhostPoints, iter, jIndex) b_acc += bΔ A_acc += AΔ end @@ -645,13 +771,14 @@ using LinearAlgebra function ApplyMDBCBeforeHalf!(SimMetaData::SimulationMetaData{D,T,S,K,SimpleMDBC,L}, SimKernel, SimConstants, SimParticles, ParticleRanges, CellDict, Position, Density, - GhostPoints, GhostNormals, ParticleType + GhostPoints, GhostNormals, ParticleType; + ParticleOrder = nothing ) where {D,T,S<:ShiftingMode,K<:KernelOutputMode,L<:LogMode} @no_escape begin 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ᵧ) + NeighborLoopMDBC!(SimKernel, SimMetaData, SimConstants, ParticleRanges, CellDict, Position, Density, GhostPoints, GhostNormals, ParticleType, bᵧ, Aᵧ; ParticleOrder = ParticleOrder) ApplyMDBCCorrection(SimConstants, SimParticles, bᵧ, Aᵧ) end @@ -709,7 +836,7 @@ using LinearAlgebra SimMetaData::SimulationMetaData{Dimensions, FloatType, SMode, KMode, BMode, LMode}, SimConstants, SimParticles, FullStencil, ParticleRanges, UniqueCells, CellDict, - SortingScratchSpace, + SortingScratchSpace, ParticleOrder, CellOffsets, NeighborCellLists, dρdtI, Velocityₙ⁺, Positionₙ⁺, ρₙ⁺, ∇Cᵢ, ∇◌rᵢ, MotionDefinition::Union{ @@ -758,7 +885,11 @@ 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⁻¹, SortingScratchSpace, ParticleRanges, UniqueCells, CellDict) + if SimMetaData.UseSortingNeighborList + @timeit SimMetaData.HourGlass "01a Actual Calculate IndexCounter" SimMetaData.IndexCounter = UpdateNeighbors!(SimParticles, SimKernel.H⁻¹, SortingScratchSpace, ParticleRanges, UniqueCells, CellDict) + else + @timeit SimMetaData.HourGlass "01a Actual Calculate IndexCounter" SimMetaData.IndexCounter = UpdateNeighborsNoSort!(SimParticles, SimKernel.H⁻¹, ParticleRanges, UniqueCells, CellDict, ParticleOrder, CellOffsets) + end SimMetaData.Δx = zero(eltype(dρdtI)) UniqueCellsView = view(UniqueCells, 1:SimMetaData.IndexCounter) BuildNeighborCellLists!(NeighborCellLists, FullStencil, UniqueCellsView, ParticleRanges, CellDict) @@ -768,10 +899,11 @@ using LinearAlgebra @timeit SimMetaData.HourGlass "Motion" ProgressMotion(SimParticles, dt₂, MotionDefinition, SimMetaData) @timeit SimMetaData.HourGlass "02 Pressure" Pressure!(SimParticles.Pressure,SimParticles.Density,SimConstants) + ParticleOrderView = SimMetaData.UseSortingNeighborList ? nothing : ParticleOrder 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) + @timeit SimMetaData.HourGlass "03 Apply MDBC before Half TimeStep" ApplyMDBCBeforeHalf!(SimMetaData, SimKernel, SimConstants, SimParticles, ParticleRanges, CellDict, Position, Density, GhostPoints, GhostNormals, ParticleType; ParticleOrder = ParticleOrderView) end if SimMetaData isa SimulationMetaData{Dimensions, FloatType, SMode, NoKernelOutput, BMode, LMode} where {SMode, BMode, LMode} @@ -779,12 +911,14 @@ using LinearAlgebra SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData, SimConstants, SimParticles, ParticleRanges, CellDict, NeighborCellLists, dρdtI, Acceleration, ∇Cᵢ, ∇◌rᵢ, + ParticleOrder = ParticleOrderView, ) else @timeit SimMetaData.HourGlass "04 First NeighborLoop" NeighborLoopPerParticle!( SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData, SimConstants, SimParticles, ParticleRanges, CellDict, NeighborCellLists, dρdtI, Acceleration, ∇Cᵢ, ∇◌rᵢ, + ParticleOrder = ParticleOrderView, ) end @@ -803,6 +937,7 @@ using LinearAlgebra SimConstants, SimParticles, ParticleRanges, CellDict, NeighborCellLists, dρdtI, Acceleration, ∇Cᵢ, ∇◌rᵢ, max_visc, min_dt_force, + ParticleOrder = ParticleOrderView, Position = Positionₙ⁺, Density = ρₙ⁺, Velocity = Velocityₙ⁺, @@ -813,6 +948,7 @@ using LinearAlgebra SimConstants, SimParticles, ParticleRanges, CellDict, NeighborCellLists, dρdtI, Acceleration, ∇Cᵢ, ∇◌rᵢ, max_visc, min_dt_force, + ParticleOrder = ParticleOrderView, Position = Positionₙ⁺, Density = ρₙ⁺, Velocity = Velocityₙ⁺, @@ -862,6 +998,8 @@ using LinearAlgebra CellDict = Dict{CartesianIndex{Dimensions}, Int}() FullStencil = ConstructStencil(Val(Dimensions)) NeighborCellLists = [Int[] for _ in 1:length(UniqueCells)] + ParticleOrder = zeros(Int, NumberOfPoints) + CellOffsets = zeros(Int, length(ParticleRanges)) _, SortingScratchSpace = Base.Sort.make_scratch(nothing, eltype(SimParticles), NumberOfPoints) output = SetupVTKOutput(SimMetaData, SimParticles, SimKernel, Dimensions) @@ -901,6 +1039,7 @@ using LinearAlgebra SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData, SimConstants, SimParticles, FullStencil, ParticleRanges, UniqueCells, CellDict, SortingScratchSpace, + ParticleOrder, CellOffsets, NeighborCellLists, dρdtI, Velocityₙ⁺, Positionₙ⁺, ρₙ⁺, ∇Cᵢ, ∇◌rᵢ, MotionDefinition, ) diff --git a/src/SPHExample.jl b/src/SPHExample.jl index 28ef3ae7..8526c77d 100644 --- a/src/SPHExample.jl +++ b/src/SPHExample.jl @@ -57,7 +57,8 @@ module SPHExample export SimulationConstants using .SPHNeighborList - export ConstructStencil, ExtractCells!, UpdateNeighbors!, BuildNeighborCellLists!, ComputeCellParticleCounts, ComputeCellNeighborCounts + export ConstructStencil, ExtractCells!, UpdateNeighbors!, UpdateNeighborsNoSort!, + BuildNeighborCellLists!, ComputeCellParticleCounts, ComputeCellNeighborCounts using .SPHCellList export NeighborLoop!, ComputeInteractions!, RunSimulation diff --git a/src/SPHNeighborList.jl b/src/SPHNeighborList.jl index 48d245ea..599ad35c 100644 --- a/src/SPHNeighborList.jl +++ b/src/SPHNeighborList.jl @@ -1,6 +1,7 @@ module SPHNeighborList -export ConstructStencil, ExtractCells!, UpdateNeighbors!, BuildNeighborCellLists!, ComputeCellParticleCounts, ComputeCellNeighborCounts, UpdateΔx! +export ConstructStencil, ExtractCells!, UpdateNeighbors!, UpdateNeighborsNoSort!, + BuildNeighborCellLists!, ComputeCellParticleCounts, ComputeCellNeighborCounts, UpdateΔx! using StaticArrays @@ -102,6 +103,53 @@ function UpdateNeighbors!(Particles, InverseCutOff, SortingScratchSpace, return IndexCounter end +""" +Updates the neighbor list without sorting particle storage. + +This builds a per-cell particle ordering buffer so cell ranges can be iterated +without reordering the particle arrays. +""" +function UpdateNeighborsNoSort!(Particles, InverseCutOff, ParticleRanges, + UniqueCells, CellDict, ParticleOrder, CellOffsets) + ExtractCells!(Particles, InverseCutOff) + + Cells = @views Particles.Cells + @. ParticleRanges = zero(eltype(ParticleRanges)) + ParticleRanges[1] = 1 + IndexCounter = 1 + empty!(CellDict) + + @inbounds for Index in eachindex(Cells) + Cell = Cells[Index] + CellIndex = get(CellDict, Cell, 0) + if CellIndex == 0 + IndexCounter += 1 + CellIndex = IndexCounter + CellDict[Cell] = CellIndex + UniqueCells[CellIndex] = Cell + end + ParticleRanges[CellIndex] += 1 + end + + RunningIndex = 1 + @inbounds for CellIndex in 2:IndexCounter + Count = ParticleRanges[CellIndex] + ParticleRanges[CellIndex] = RunningIndex + CellOffsets[CellIndex] = RunningIndex + RunningIndex += Count + end + ParticleRanges[IndexCounter + 1] = RunningIndex + + @inbounds for Index in eachindex(Cells) + CellIndex = CellDict[Cells[Index]] + TargetIndex = CellOffsets[CellIndex] + ParticleOrder[TargetIndex] = Index + CellOffsets[CellIndex] = TargetIndex + 1 + end + + return IndexCounter +end + function ComputeCellParticleCounts(ParticleRanges, CellCount) Counts = Vector{Int}(undef, CellCount) @inbounds for Index in 1:CellCount diff --git a/src/SimulationMetaDataConfiguration.jl b/src/SimulationMetaDataConfiguration.jl index c0ce65df..56ced401 100644 --- a/src/SimulationMetaDataConfiguration.jl +++ b/src/SimulationMetaDataConfiguration.jl @@ -48,6 +48,7 @@ struct StoreLog <: LogMode end ExportSingleVTKHDF::Bool = true ExportGridCells::Bool = false ExportGridCellParticleCounts::Bool = false + UseSortingNeighborList::Bool = true OutputVariables::Vector{String} = [ "Density", "Pressure", From 6059e3c59135e3c5f8526949a582b5ea20ce3838 Mon Sep 17 00:00:00 2001 From: AhmedSalih3d <36305327+AhmedSalih3d@users.noreply.github.com> Date: Wed, 14 Jan 2026 22:33:51 +0100 Subject: [PATCH 2/8] Default to no-sort neighbor list --- README.md | 6 +++--- src/SPHNeighborList.jl | 4 ++-- src/SimulationMetaDataConfiguration.jl | 2 +- 3 files changed, 6 insertions(+), 6 deletions(-) diff --git a/README.md b/README.md index fb3cc80c..aca23a59 100644 --- a/README.md +++ b/README.md @@ -88,9 +88,9 @@ To color exported cell grids by particle counts, set `ExportGridCellParticleCounts=true` in `SimulationMetaData`. This also adds a `ParticleNeighborsPerCell` array that includes each cell's particle count minus one plus the particles in its neighbor stencil. -If you want to avoid re-sorting particles each neighbor rebuild, set -`UseSortingNeighborList=false` to build a per-cell particle ordering buffer -without reordering particle storage. +To avoid re-sorting particles each neighbor rebuild, keep +`UseSortingNeighborList=false` (the default) to build a per-cell particle +ordering buffer without reordering particle storage. ## Help diff --git a/src/SPHNeighborList.jl b/src/SPHNeighborList.jl index 599ad35c..aa4617ac 100644 --- a/src/SPHNeighborList.jl +++ b/src/SPHNeighborList.jl @@ -82,7 +82,7 @@ function UpdateNeighbors!(Particles, InverseCutOff, SortingScratchSpace, sort!(Particles, by = p -> p.Cells; scratch=SortingScratchSpace) Cells = @views Particles.Cells - @. ParticleRanges = zero(eltype(ParticleRanges)) + fill!(ParticleRanges, zero(eltype(ParticleRanges))) ParticleRanges[1] = 1 IndexCounter = 2 ParticleRanges[IndexCounter] = 1 @@ -114,7 +114,7 @@ function UpdateNeighborsNoSort!(Particles, InverseCutOff, ParticleRanges, ExtractCells!(Particles, InverseCutOff) Cells = @views Particles.Cells - @. ParticleRanges = zero(eltype(ParticleRanges)) + fill!(ParticleRanges, zero(eltype(ParticleRanges))) ParticleRanges[1] = 1 IndexCounter = 1 empty!(CellDict) diff --git a/src/SimulationMetaDataConfiguration.jl b/src/SimulationMetaDataConfiguration.jl index 56ced401..97d21f41 100644 --- a/src/SimulationMetaDataConfiguration.jl +++ b/src/SimulationMetaDataConfiguration.jl @@ -48,7 +48,7 @@ struct StoreLog <: LogMode end ExportSingleVTKHDF::Bool = true ExportGridCells::Bool = false ExportGridCellParticleCounts::Bool = false - UseSortingNeighborList::Bool = true + UseSortingNeighborList::Bool = false OutputVariables::Vector{String} = [ "Density", "Pressure", From bfcd68010fee361a2522baffbafc61f0f2afcae1 Mon Sep 17 00:00:00 2001 From: AhmedSalih3d <36305327+AhmedSalih3d@users.noreply.github.com> Date: Wed, 14 Jan 2026 22:37:56 +0100 Subject: [PATCH 3/8] Simplify to single neighbor list path --- README.md | 5 +- src/SPHCellList.jl | 267 ++++++------------------- src/SPHExample.jl | 2 +- src/SPHNeighborList.jl | 46 +---- src/SimulationMetaDataConfiguration.jl | 1 - 5 files changed, 67 insertions(+), 254 deletions(-) diff --git a/README.md b/README.md index aca23a59..c23bc265 100644 --- a/README.md +++ b/README.md @@ -88,9 +88,8 @@ To color exported cell grids by particle counts, set `ExportGridCellParticleCounts=true` in `SimulationMetaData`. This also adds a `ParticleNeighborsPerCell` array that includes each cell's particle count minus one plus the particles in its neighbor stencil. -To avoid re-sorting particles each neighbor rebuild, keep -`UseSortingNeighborList=false` (the default) to build a per-cell particle -ordering buffer without reordering particle storage. +Neighbor rebuilds use a per-cell ordering buffer without reordering particle +storage. ## Help diff --git a/src/SPHCellList.jl b/src/SPHCellList.jl index feedc1b8..2076f2cc 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!, UpdateNeighborsNoSort!, UpdateΔx! +using ..SPHNeighborList: BuildNeighborCellLists!, ComputeCellNeighborCounts, ComputeCellParticleCounts, ConstructStencil, ExtractCells!, MapFloor, UpdateNeighbors!, UpdateΔx! using StaticArrays using StructArrays: StructArray, foreachfield @@ -32,14 +32,6 @@ using Base.Threads using LinearAlgebra using Bumper - @inline function ResolveParticleIndex(::Nothing, Index) - return Index - end - - @inline function ResolveParticleIndex(ParticleOrder, Index) - return @inbounds ParticleOrder[Index] - end - function NeighborLoopPerParticle!(SimDensityDiffusion::SDD, SimViscosity::SV, SimKernel, SimMetaData::SimulationMetaData{D,T,NoShifting,NoKernelOutput,B,L}, SimConstants, SimParticles, ParticleRanges, @@ -51,7 +43,7 @@ using LinearAlgebra Density = SimParticles.Density, Pressure = SimParticles.Pressure, Velocity = SimParticles.Velocity, - ParticleOrder = nothing) where {D,T, + ParticleOrder) where {D,T, B<:MDBCMode,L<:LogMode, SDD<:SPHDensityDiffusion, SV<:SPHViscosity} @@ -65,55 +57,27 @@ using LinearAlgebra SameCellEnd = ParticleRanges[CellListIndex + 1] - 1 NeighborCellIndices = NeighborCellLists[CellListIndex] - if ParticleOrder === nothing - @inbounds for j in SameCellStart:(i - 1) + @inbounds for j in SameCellStart:SameCellEnd + jIndex = ParticleOrder[j] + if jIndex != i dρdt_acc, acc_acc = ComputeInteractionsPerParticleNoKernel!( SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData, SimConstants, SimParticles, Position, Density, Pressure, - Velocity, MotionLimiter, dρdt_acc, acc_acc, i, j, + Velocity, MotionLimiter, dρdt_acc, acc_acc, i, jIndex, ) end - @inbounds for j in (i + 1):SameCellEnd + end + for NeighborIdx in NeighborCellIndices + StartIndex_ = ParticleRanges[NeighborIdx] + EndIndex_ = ParticleRanges[NeighborIdx + 1] - 1 + @inbounds for j in StartIndex_:EndIndex_ + jIndex = ParticleOrder[j] dρdt_acc, acc_acc = ComputeInteractionsPerParticleNoKernel!( SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData, SimConstants, SimParticles, Position, Density, Pressure, - Velocity, MotionLimiter, dρdt_acc, acc_acc, i, j, + Velocity, MotionLimiter, dρdt_acc, acc_acc, i, jIndex, ) end - for NeighborIdx in NeighborCellIndices - StartIndex_ = ParticleRanges[NeighborIdx] - EndIndex_ = ParticleRanges[NeighborIdx + 1] - 1 - @inbounds for j in StartIndex_:EndIndex_ - dρdt_acc, acc_acc = ComputeInteractionsPerParticleNoKernel!( - SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData, - SimConstants, SimParticles, Position, Density, Pressure, - Velocity, MotionLimiter, dρdt_acc, acc_acc, i, j, - ) - end - end - else - @inbounds for j in SameCellStart:SameCellEnd - jIndex = ResolveParticleIndex(ParticleOrder, j) - if jIndex != i - dρdt_acc, acc_acc = ComputeInteractionsPerParticleNoKernel!( - SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData, - SimConstants, SimParticles, Position, Density, Pressure, - Velocity, MotionLimiter, dρdt_acc, acc_acc, i, jIndex, - ) - end - end - for NeighborIdx in NeighborCellIndices - StartIndex_ = ParticleRanges[NeighborIdx] - EndIndex_ = ParticleRanges[NeighborIdx + 1] - 1 - @inbounds for j in StartIndex_:EndIndex_ - jIndex = ResolveParticleIndex(ParticleOrder, j) - dρdt_acc, acc_acc = ComputeInteractionsPerParticleNoKernel!( - SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData, - SimConstants, SimParticles, Position, Density, Pressure, - Velocity, MotionLimiter, dρdt_acc, acc_acc, i, jIndex, - ) - end - end end dρdtI[i] = dρdt_acc @@ -135,7 +99,7 @@ using LinearAlgebra Density = SimParticles.Density, Pressure = SimParticles.Pressure, Velocity = SimParticles.Velocity, - ParticleOrder = nothing) where {D,T, + ParticleOrder) where {D,T, K<:KernelOutputMode, B<:MDBCMode,L<:LogMode, SDD<:SPHDensityDiffusion, @@ -152,65 +116,31 @@ using LinearAlgebra SameCellEnd = ParticleRanges[CellListIndex + 1] - 1 NeighborCellIndices = NeighborCellLists[CellListIndex] - if ParticleOrder === nothing - @inbounds for j in SameCellStart:(i - 1) + @inbounds for j in SameCellStart:SameCellEnd + jIndex = ParticleOrder[j] + if jIndex != i dρdt_acc, acc_acc, kernel_acc, kernel_grad_acc = ComputeInteractionsPerParticle!( SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData, SimConstants, SimParticles, Position, Density, Pressure, Velocity, MotionLimiter, dρdt_acc, acc_acc, kernel_acc, - kernel_grad_acc, i, j, + kernel_grad_acc, i, jIndex, ) end - @inbounds for j in (i + 1):SameCellEnd + end + for NeighborIdx in NeighborCellIndices + StartIndex_ = ParticleRanges[NeighborIdx] + EndIndex_ = ParticleRanges[NeighborIdx + 1] - 1 + @inbounds for j in StartIndex_:EndIndex_ + jIndex = ParticleOrder[j] dρdt_acc, acc_acc, kernel_acc, kernel_grad_acc = ComputeInteractionsPerParticle!( SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData, SimConstants, SimParticles, Position, Density, Pressure, Velocity, MotionLimiter, dρdt_acc, acc_acc, kernel_acc, - kernel_grad_acc, i, j, + kernel_grad_acc, i, jIndex, ) end - for NeighborIdx in NeighborCellIndices - StartIndex_ = ParticleRanges[NeighborIdx] - EndIndex_ = ParticleRanges[NeighborIdx + 1] - 1 - @inbounds for j in StartIndex_:EndIndex_ - dρdt_acc, acc_acc, kernel_acc, kernel_grad_acc = - ComputeInteractionsPerParticle!( - SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData, - SimConstants, SimParticles, Position, Density, Pressure, - Velocity, MotionLimiter, dρdt_acc, acc_acc, kernel_acc, - kernel_grad_acc, i, j, - ) - end - end - else - @inbounds for j in SameCellStart:SameCellEnd - jIndex = ResolveParticleIndex(ParticleOrder, j) - if jIndex != i - dρdt_acc, acc_acc, kernel_acc, kernel_grad_acc = - ComputeInteractionsPerParticle!( - SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData, - SimConstants, SimParticles, Position, Density, Pressure, - Velocity, MotionLimiter, dρdt_acc, acc_acc, kernel_acc, - kernel_grad_acc, i, jIndex, - ) - end - end - for NeighborIdx in NeighborCellIndices - StartIndex_ = ParticleRanges[NeighborIdx] - EndIndex_ = ParticleRanges[NeighborIdx + 1] - 1 - @inbounds for j in StartIndex_:EndIndex_ - jIndex = ResolveParticleIndex(ParticleOrder, j) - dρdt_acc, acc_acc, kernel_acc, kernel_grad_acc = - ComputeInteractionsPerParticle!( - SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData, - SimConstants, SimParticles, Position, Density, Pressure, - Velocity, MotionLimiter, dρdt_acc, acc_acc, kernel_acc, - kernel_grad_acc, i, jIndex, - ) - end - end end dρdtI[i] = dρdt_acc @@ -234,7 +164,7 @@ using LinearAlgebra Density = SimParticles.Density, Pressure = SimParticles.Pressure, Velocity = SimParticles.Velocity, - ParticleOrder = nothing) where {D,T, + ParticleOrder) where {D,T, S<:ShiftingMode,B<:MDBCMode, L<:LogMode,SDD<:SPHDensityDiffusion, SV<:SPHViscosity} @@ -250,65 +180,31 @@ using LinearAlgebra SameCellEnd = ParticleRanges[CellListIndex + 1] - 1 NeighborCellIndices = NeighborCellLists[CellListIndex] - if ParticleOrder === nothing - @inbounds for j in SameCellStart:(i - 1) + @inbounds for j in SameCellStart:SameCellEnd + jIndex = ParticleOrder[j] + if jIndex != i dρdt_acc, acc_acc, shift_c_acc, shift_r_acc = ComputeInteractionsPerParticleNoKernel!( SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData, SimConstants, SimParticles, Position, Density, Pressure, Velocity, MotionLimiter, dρdt_acc, acc_acc, shift_c_acc, - shift_r_acc, i, j, + shift_r_acc, i, jIndex, ) end - @inbounds for j in (i + 1):SameCellEnd + end + for NeighborIdx in NeighborCellIndices + StartIndex_ = ParticleRanges[NeighborIdx] + EndIndex_ = ParticleRanges[NeighborIdx + 1] - 1 + @inbounds for j in StartIndex_:EndIndex_ + jIndex = ParticleOrder[j] dρdt_acc, acc_acc, shift_c_acc, shift_r_acc = ComputeInteractionsPerParticleNoKernel!( SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData, SimConstants, SimParticles, Position, Density, Pressure, Velocity, MotionLimiter, dρdt_acc, acc_acc, shift_c_acc, - shift_r_acc, i, j, + shift_r_acc, i, jIndex, ) end - for NeighborIdx in NeighborCellIndices - StartIndex_ = ParticleRanges[NeighborIdx] - EndIndex_ = ParticleRanges[NeighborIdx + 1] - 1 - @inbounds for j in StartIndex_:EndIndex_ - dρdt_acc, acc_acc, shift_c_acc, shift_r_acc = - ComputeInteractionsPerParticleNoKernel!( - SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData, - SimConstants, SimParticles, Position, Density, Pressure, - Velocity, MotionLimiter, dρdt_acc, acc_acc, shift_c_acc, - shift_r_acc, i, j, - ) - end - end - else - @inbounds for j in SameCellStart:SameCellEnd - jIndex = ResolveParticleIndex(ParticleOrder, j) - if jIndex != i - dρdt_acc, acc_acc, shift_c_acc, shift_r_acc = - ComputeInteractionsPerParticleNoKernel!( - SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData, - SimConstants, SimParticles, Position, Density, Pressure, - Velocity, MotionLimiter, dρdt_acc, acc_acc, shift_c_acc, - shift_r_acc, i, jIndex, - ) - end - end - for NeighborIdx in NeighborCellIndices - StartIndex_ = ParticleRanges[NeighborIdx] - EndIndex_ = ParticleRanges[NeighborIdx + 1] - 1 - @inbounds for j in StartIndex_:EndIndex_ - jIndex = ResolveParticleIndex(ParticleOrder, j) - dρdt_acc, acc_acc, shift_c_acc, shift_r_acc = - ComputeInteractionsPerParticleNoKernel!( - SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData, - SimConstants, SimParticles, Position, Density, Pressure, - Velocity, MotionLimiter, dρdt_acc, acc_acc, shift_c_acc, - shift_r_acc, i, jIndex, - ) - end - end end dρdtI[i] = dρdt_acc @@ -332,7 +228,7 @@ using LinearAlgebra Density = SimParticles.Density, Pressure = SimParticles.Pressure, Velocity = SimParticles.Velocity, - ParticleOrder = nothing) where {D,T, + ParticleOrder) where {D,T, S<:ShiftingMode, K<:KernelOutputMode, B<:MDBCMode,L<:LogMode, @@ -352,65 +248,31 @@ using LinearAlgebra SameCellEnd = ParticleRanges[CellListIndex + 1] - 1 NeighborCellIndices = NeighborCellLists[CellListIndex] - if ParticleOrder === nothing - @inbounds for j in SameCellStart:(i - 1) + @inbounds for j in SameCellStart:SameCellEnd + jIndex = ParticleOrder[j] + if jIndex != i dρdt_acc, acc_acc, kernel_acc, kernel_grad_acc, shift_c_acc, shift_r_acc = ComputeInteractionsPerParticle!( SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData, SimConstants, SimParticles, Position, Density, Pressure, Velocity, MotionLimiter, dρdt_acc, acc_acc, kernel_acc, - kernel_grad_acc, shift_c_acc, shift_r_acc, i, j, + kernel_grad_acc, shift_c_acc, shift_r_acc, i, jIndex, ) end - @inbounds for j in (i + 1):SameCellEnd + end + for NeighborIdx in NeighborCellIndices + StartIndex_ = ParticleRanges[NeighborIdx] + EndIndex_ = ParticleRanges[NeighborIdx + 1] - 1 + @inbounds for j in StartIndex_:EndIndex_ + jIndex = ParticleOrder[j] dρdt_acc, acc_acc, kernel_acc, kernel_grad_acc, shift_c_acc, shift_r_acc = ComputeInteractionsPerParticle!( SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData, SimConstants, SimParticles, Position, Density, Pressure, Velocity, MotionLimiter, dρdt_acc, acc_acc, kernel_acc, - kernel_grad_acc, shift_c_acc, shift_r_acc, i, j, + kernel_grad_acc, shift_c_acc, shift_r_acc, i, jIndex, ) end - for NeighborIdx in NeighborCellIndices - StartIndex_ = ParticleRanges[NeighborIdx] - EndIndex_ = ParticleRanges[NeighborIdx + 1] - 1 - @inbounds for j in StartIndex_:EndIndex_ - dρdt_acc, acc_acc, kernel_acc, kernel_grad_acc, shift_c_acc, - shift_r_acc = ComputeInteractionsPerParticle!( - SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData, - SimConstants, SimParticles, Position, Density, Pressure, - Velocity, MotionLimiter, dρdt_acc, acc_acc, kernel_acc, - kernel_grad_acc, shift_c_acc, shift_r_acc, i, j, - ) - end - end - else - @inbounds for j in SameCellStart:SameCellEnd - jIndex = ResolveParticleIndex(ParticleOrder, j) - if jIndex != i - dρdt_acc, acc_acc, kernel_acc, kernel_grad_acc, shift_c_acc, - shift_r_acc = ComputeInteractionsPerParticle!( - SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData, - SimConstants, SimParticles, Position, Density, Pressure, - Velocity, MotionLimiter, dρdt_acc, acc_acc, kernel_acc, - kernel_grad_acc, shift_c_acc, shift_r_acc, i, jIndex, - ) - end - end - for NeighborIdx in NeighborCellIndices - StartIndex_ = ParticleRanges[NeighborIdx] - EndIndex_ = ParticleRanges[NeighborIdx + 1] - 1 - @inbounds for j in StartIndex_:EndIndex_ - jIndex = ResolveParticleIndex(ParticleOrder, j) - dρdt_acc, acc_acc, kernel_acc, kernel_grad_acc, shift_c_acc, - shift_r_acc = ComputeInteractionsPerParticle!( - SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData, - SimConstants, SimParticles, Position, Density, Pressure, - Velocity, MotionLimiter, dρdt_acc, acc_acc, kernel_acc, - kernel_grad_acc, shift_c_acc, shift_r_acc, i, jIndex, - ) - end - end end dρdtI[i] = dρdt_acc @@ -432,7 +294,7 @@ using LinearAlgebra SimConstants, ParticleRanges, CellDict, Position, Density, GhostPoints, GhostNormals, ParticleType, bᵧ, Aᵧ; - ParticleOrder = nothing) where {Dimensions, FloatType, SMode, KMode, BMode, LMode} + ParticleOrder) where {Dimensions, FloatType, SMode, KMode, BMode, LMode} FullStencil = ConstructStencil(Val(Dimensions)) @@ -457,9 +319,9 @@ using LinearAlgebra StartIndex_ = ParticleRanges[NeighborIdx] EndIndex_ = ParticleRanges[NeighborIdx + 1] - 1 - for j in StartIndex_:EndIndex_ + @inbounds for j in StartIndex_:EndIndex_ # change ComputeInteractions to take & return contributions, e.g.: - jIndex = ResolveParticleIndex(ParticleOrder, j) + jIndex = ParticleOrder[j] bΔ, AΔ = ComputeInteractionsMDBC!(SimKernel, SimMetaData, SimConstants, Position, Density, ParticleType, GhostPoints, iter, jIndex) @@ -772,7 +634,7 @@ using LinearAlgebra SimKernel, SimConstants, SimParticles, ParticleRanges, CellDict, Position, Density, GhostPoints, GhostNormals, ParticleType; - ParticleOrder = nothing + ParticleOrder ) where {D,T,S<:ShiftingMode,K<:KernelOutputMode,L<:LogMode} @no_escape begin DimensionsPlus = D + 1 @@ -836,7 +698,7 @@ using LinearAlgebra SimMetaData::SimulationMetaData{Dimensions, FloatType, SMode, KMode, BMode, LMode}, SimConstants, SimParticles, FullStencil, ParticleRanges, UniqueCells, CellDict, - SortingScratchSpace, ParticleOrder, CellOffsets, + ParticleOrder, CellOffsets, NeighborCellLists, dρdtI, Velocityₙ⁺, Positionₙ⁺, ρₙ⁺, ∇Cᵢ, ∇◌rᵢ, MotionDefinition::Union{ @@ -885,11 +747,7 @@ 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 - if SimMetaData.UseSortingNeighborList - @timeit SimMetaData.HourGlass "01a Actual Calculate IndexCounter" SimMetaData.IndexCounter = UpdateNeighbors!(SimParticles, SimKernel.H⁻¹, SortingScratchSpace, ParticleRanges, UniqueCells, CellDict) - else - @timeit SimMetaData.HourGlass "01a Actual Calculate IndexCounter" SimMetaData.IndexCounter = UpdateNeighborsNoSort!(SimParticles, SimKernel.H⁻¹, ParticleRanges, UniqueCells, CellDict, ParticleOrder, CellOffsets) - end + @timeit SimMetaData.HourGlass "01a Actual Calculate IndexCounter" SimMetaData.IndexCounter = UpdateNeighbors!(SimParticles, SimKernel.H⁻¹, ParticleRanges, UniqueCells, CellDict, ParticleOrder, CellOffsets) SimMetaData.Δx = zero(eltype(dρdtI)) UniqueCellsView = view(UniqueCells, 1:SimMetaData.IndexCounter) BuildNeighborCellLists!(NeighborCellLists, FullStencil, UniqueCellsView, ParticleRanges, CellDict) @@ -899,11 +757,10 @@ using LinearAlgebra @timeit SimMetaData.HourGlass "Motion" ProgressMotion(SimParticles, dt₂, MotionDefinition, SimMetaData) @timeit SimMetaData.HourGlass "02 Pressure" Pressure!(SimParticles.Pressure,SimParticles.Density,SimConstants) - ParticleOrderView = SimMetaData.UseSortingNeighborList ? nothing : ParticleOrder 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 = ParticleOrderView) + @timeit SimMetaData.HourGlass "03 Apply MDBC before Half TimeStep" ApplyMDBCBeforeHalf!(SimMetaData, SimKernel, SimConstants, SimParticles, ParticleRanges, CellDict, Position, Density, GhostPoints, GhostNormals, ParticleType; ParticleOrder = ParticleOrder) end if SimMetaData isa SimulationMetaData{Dimensions, FloatType, SMode, NoKernelOutput, BMode, LMode} where {SMode, BMode, LMode} @@ -911,14 +768,14 @@ using LinearAlgebra SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData, SimConstants, SimParticles, ParticleRanges, CellDict, NeighborCellLists, dρdtI, Acceleration, ∇Cᵢ, ∇◌rᵢ, - ParticleOrder = ParticleOrderView, + ParticleOrder = ParticleOrder, ) else @timeit SimMetaData.HourGlass "04 First NeighborLoop" NeighborLoopPerParticle!( SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData, SimConstants, SimParticles, ParticleRanges, CellDict, NeighborCellLists, dρdtI, Acceleration, ∇Cᵢ, ∇◌rᵢ, - ParticleOrder = ParticleOrderView, + ParticleOrder = ParticleOrder, ) end @@ -937,7 +794,7 @@ using LinearAlgebra SimConstants, SimParticles, ParticleRanges, CellDict, NeighborCellLists, dρdtI, Acceleration, ∇Cᵢ, ∇◌rᵢ, max_visc, min_dt_force, - ParticleOrder = ParticleOrderView, + ParticleOrder = ParticleOrder, Position = Positionₙ⁺, Density = ρₙ⁺, Velocity = Velocityₙ⁺, @@ -948,7 +805,7 @@ using LinearAlgebra SimConstants, SimParticles, ParticleRanges, CellDict, NeighborCellLists, dρdtI, Acceleration, ∇Cᵢ, ∇◌rᵢ, max_visc, min_dt_force, - ParticleOrder = ParticleOrderView, + ParticleOrder = ParticleOrder, Position = Positionₙ⁺, Density = ρₙ⁺, Velocity = Velocityₙ⁺, @@ -1000,7 +857,6 @@ using LinearAlgebra NeighborCellLists = [Int[] for _ in 1:length(UniqueCells)] ParticleOrder = zeros(Int, NumberOfPoints) CellOffsets = zeros(Int, length(ParticleRanges)) - _, SortingScratchSpace = Base.Sort.make_scratch(nothing, eltype(SimParticles), NumberOfPoints) output = SetupVTKOutput(SimMetaData, SimParticles, SimKernel, Dimensions) @@ -1038,8 +894,7 @@ using LinearAlgebra @timeit SimMetaData.HourGlass "00 SimulationLoop" SimulationLoop( SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData, SimConstants, SimParticles, FullStencil, ParticleRanges, - UniqueCells, CellDict, SortingScratchSpace, - ParticleOrder, CellOffsets, + UniqueCells, CellDict, ParticleOrder, CellOffsets, NeighborCellLists, dρdtI, Velocityₙ⁺, Positionₙ⁺, ρₙ⁺, ∇Cᵢ, ∇◌rᵢ, MotionDefinition, ) diff --git a/src/SPHExample.jl b/src/SPHExample.jl index 8526c77d..a56895ea 100644 --- a/src/SPHExample.jl +++ b/src/SPHExample.jl @@ -57,7 +57,7 @@ module SPHExample export SimulationConstants using .SPHNeighborList - export ConstructStencil, ExtractCells!, UpdateNeighbors!, UpdateNeighborsNoSort!, + export ConstructStencil, ExtractCells!, UpdateNeighbors!, BuildNeighborCellLists!, ComputeCellParticleCounts, ComputeCellNeighborCounts using .SPHCellList diff --git a/src/SPHNeighborList.jl b/src/SPHNeighborList.jl index aa4617ac..39c4c2a0 100644 --- a/src/SPHNeighborList.jl +++ b/src/SPHNeighborList.jl @@ -1,6 +1,6 @@ module SPHNeighborList -export ConstructStencil, ExtractCells!, UpdateNeighbors!, UpdateNeighborsNoSort!, +export ConstructStencil, ExtractCells!, UpdateNeighbors!, BuildNeighborCellLists!, ComputeCellParticleCounts, ComputeCellNeighborCounts, UpdateΔx! using StaticArrays @@ -63,54 +63,14 @@ end return nothing end -""" -Updates the neighbor list and sorts particles by their cell indices. - -# Arguments -- `Particles`: The particles whose neighbors are to be updated. -- `CutOff`: The cutoff value used for cell extraction. -- `SortingScratchSpace`: Scratch space for sorting. -- `ParticleRanges`: Array to store the ranges of particles in each cell. -- `UniqueCells`: Array to store the unique cells. - -# Returns -- `IndexCounter`: The number of unique cells identified. -""" -function UpdateNeighbors!(Particles, InverseCutOff, SortingScratchSpace, - ParticleRanges, UniqueCells, CellDict) - ExtractCells!(Particles, InverseCutOff) - - sort!(Particles, by = p -> p.Cells; scratch=SortingScratchSpace) - Cells = @views Particles.Cells - fill!(ParticleRanges, zero(eltype(ParticleRanges))) - ParticleRanges[1] = 1 - IndexCounter = 2 - ParticleRanges[IndexCounter] = 1 - UniqueCells[IndexCounter] = Cells[1] - empty!(CellDict) - CellDict[Cells[1]] = IndexCounter - - @inbounds @simd ivdep for Index in eachindex(Cells)[2:end] - if Cells[Index] != Cells[Index - 1] # Equivalent to diff(Cells) != 0 - IndexCounter += 1 - ParticleRanges[IndexCounter] = Index - UniqueCells[IndexCounter] = Cells[Index] - CellDict[Cells[Index]] = IndexCounter - end - end - ParticleRanges[IndexCounter + 1] = length(ParticleRanges) - - return IndexCounter -end - """ Updates the neighbor list without sorting particle storage. This builds a per-cell particle ordering buffer so cell ranges can be iterated without reordering the particle arrays. """ -function UpdateNeighborsNoSort!(Particles, InverseCutOff, ParticleRanges, - UniqueCells, CellDict, ParticleOrder, CellOffsets) +function UpdateNeighbors!(Particles, InverseCutOff, ParticleRanges, + UniqueCells, CellDict, ParticleOrder, CellOffsets) ExtractCells!(Particles, InverseCutOff) Cells = @views Particles.Cells diff --git a/src/SimulationMetaDataConfiguration.jl b/src/SimulationMetaDataConfiguration.jl index 97d21f41..c0ce65df 100644 --- a/src/SimulationMetaDataConfiguration.jl +++ b/src/SimulationMetaDataConfiguration.jl @@ -48,7 +48,6 @@ struct StoreLog <: LogMode end ExportSingleVTKHDF::Bool = true ExportGridCells::Bool = false ExportGridCellParticleCounts::Bool = false - UseSortingNeighborList::Bool = false OutputVariables::Vector{String} = [ "Density", "Pressure", From fccfce32643a90308fc36f4223b0476c3b1b0580 Mon Sep 17 00:00:00 2001 From: AhmedSalih3d <36305327+AhmedSalih3d@users.noreply.github.com> Date: Wed, 14 Jan 2026 22:44:00 +0100 Subject: [PATCH 4/8] Speed up neighbor cell list build --- src/SPHNeighborList.jl | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 deletions(-) diff --git a/src/SPHNeighborList.jl b/src/SPHNeighborList.jl index 39c4c2a0..924e20f9 100644 --- a/src/SPHNeighborList.jl +++ b/src/SPHNeighborList.jl @@ -23,14 +23,17 @@ function BuildNeighborCellLists!(NeighborCellLists, FullStencil, UniqueCellsView @inbounds for CellIndex in eachindex(UniqueCellsView) Neighbors = NeighborCellLists[CellIndex] empty!(Neighbors) + sizehint!(Neighbors, length(FullStencil) - 1) Cell = UniqueCellsView[CellIndex] for Offset in FullStencil NeighborCell = Cell + Offset - NeighborIndex = get(CellDict, NeighborCell, 1) - StartIndex = ParticleRanges[NeighborIndex] - EndIndex = ParticleRanges[NeighborIndex + 1] - 1 - if StartIndex <= EndIndex && NeighborIndex != CellIndex - push!(Neighbors, NeighborIndex) + 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 end end end From 78172ee5f7311e4d3e218d9dbc4c43a4747635b9 Mon Sep 17 00:00:00 2001 From: AhmedSalih3d <36305327+AhmedSalih3d@users.noreply.github.com> Date: Wed, 14 Jan 2026 22:52:34 +0100 Subject: [PATCH 5/8] Reduce neighbor list rebuild passes --- src/SPHNeighborList.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/SPHNeighborList.jl b/src/SPHNeighborList.jl index 924e20f9..b94fc07e 100644 --- a/src/SPHNeighborList.jl +++ b/src/SPHNeighborList.jl @@ -77,10 +77,10 @@ function UpdateNeighbors!(Particles, InverseCutOff, ParticleRanges, ExtractCells!(Particles, InverseCutOff) Cells = @views Particles.Cells - fill!(ParticleRanges, zero(eltype(ParticleRanges))) ParticleRanges[1] = 1 IndexCounter = 1 empty!(CellDict) + fill!(CellOffsets, zero(eltype(CellOffsets))) @inbounds for Index in eachindex(Cells) Cell = Cells[Index] @@ -91,12 +91,12 @@ function UpdateNeighbors!(Particles, InverseCutOff, ParticleRanges, CellDict[Cell] = CellIndex UniqueCells[CellIndex] = Cell end - ParticleRanges[CellIndex] += 1 + CellOffsets[CellIndex] += 1 end RunningIndex = 1 @inbounds for CellIndex in 2:IndexCounter - Count = ParticleRanges[CellIndex] + Count = CellOffsets[CellIndex] ParticleRanges[CellIndex] = RunningIndex CellOffsets[CellIndex] = RunningIndex RunningIndex += Count From 1bacdeaa7795242a73187a2fde07266f30f250db Mon Sep 17 00:00:00 2001 From: AhmedSalih3d <36305327+AhmedSalih3d@users.noreply.github.com> Date: Thu, 15 Jan 2026 20:07:56 +0100 Subject: [PATCH 6/8] Slight refactor --- src/SPHCellList.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/SPHCellList.jl b/src/SPHCellList.jl index 2076f2cc..ca55231f 100644 --- a/src/SPHCellList.jl +++ b/src/SPHCellList.jl @@ -932,11 +932,11 @@ using LinearAlgebra show(SimMetaData.HourGlass,sortby=:name) show(SimMetaData.HourGlass) - AutoOpenParaview(SimMetaData, output.variable_names) - FinalizeLog!(SimMetaData, SimLogger) + AutoOpenLogFile(SimLogger, SimMetaData) - + AutoOpenParaview(SimMetaData, output.variable_names) + break end end From 9f470bb199ad3076ff367a021541db4aea96542e Mon Sep 17 00:00:00 2001 From: AhmedSalih3d <36305327+AhmedSalih3d@users.noreply.github.com> Date: Thu, 15 Jan 2026 20:25:31 +0100 Subject: [PATCH 7/8] Using multiple dispatch for outputting grid in simulation. In shaa Allah this will make code execute faster --- src/SPHCellList.jl | 36 ++++++++++++++++++------------------ 1 file changed, 18 insertions(+), 18 deletions(-) diff --git a/src/SPHCellList.jl b/src/SPHCellList.jl index ca55231f..19ddae69 100644 --- a/src/SPHCellList.jl +++ b/src/SPHCellList.jl @@ -826,6 +826,21 @@ using LinearAlgebra return nothing 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)) + + @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) + end + end + + function DetermineOutput(::Val{false}, SimMetaData, output, ParticleRanges, NeighborCellLists, UniqueCellsView) + output.enqueue_particles(SimMetaData.OutputIterationCounter) + end ###=== function RunSimulation(;SimGeometry::Vector{Geometry{Dimensions, FloatType}}, #Don't further specify type for now @@ -905,24 +920,9 @@ using LinearAlgebra SimMetaData.OutputIterationCounter += 1 UniqueCellsView = view(UniqueCells, 1:SimMetaData.IndexCounter) - cell_particle_counts = nothing - cell_neighbor_counts = nothing - if SimMetaData.ExportGridCellParticleCounts - cell_particle_counts = ComputeCellParticleCounts( - ParticleRanges, - length(UniqueCellsView), - ) - cell_neighbor_counts = ComputeCellNeighborCounts( - ParticleRanges, - NeighborCellLists, - length(UniqueCellsView), - ) - end - - @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) - end + + @timeit SimMetaData.HourGlass "13 Determine Output" DetermineOutput(Val(SimMetaData.ExportGridCellParticleCounts), SimMetaData,output, ParticleRanges, NeighborCellLists, UniqueCellsView) + if SimMetaData.TotalTime > SimMetaData.SimulationTime From 4ec7dfb2ed20226c8d73f9a3736ac40c0c7aa55b Mon Sep 17 00:00:00 2001 From: AhmedSalih3d <36305327+AhmedSalih3d@users.noreply.github.com> Date: Thu, 15 Jan 2026 22:17:05 +0100 Subject: [PATCH 8/8] Replace dict cell lookup --- src/SPHCellList.jl | 48 +++++++++---------- src/SPHExample.jl | 3 +- src/SPHNeighborList.jl | 102 +++++++++++++++++++++++++++++++++++++---- 3 files changed, 120 insertions(+), 33 deletions(-) diff --git a/src/SPHCellList.jl b/src/SPHCellList.jl index 19ddae69..7d149e2e 100644 --- a/src/SPHCellList.jl +++ b/src/SPHCellList.jl @@ -18,7 +18,9 @@ using ..OpenExternalPrograms using ..SPHKernels using ..SPHViscosityModels using ..SPHDensityDiffusionModels -using ..SPHNeighborList: BuildNeighborCellLists!, ComputeCellNeighborCounts, ComputeCellParticleCounts, ConstructStencil, ExtractCells!, MapFloor, UpdateNeighbors!, UpdateΔx! +using ..SPHNeighborList: BuildNeighborCellLists!, CellLookupIndex, ComputeCellNeighborCounts, + ComputeCellParticleCounts, ConstructStencil, ExtractCells!, InitializeCellIndexLookup, + 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, - CellDict, NeighborCellLists, dρdtI, + CellLookup, NeighborCellLists, dρdtI, Acceleration, ∇Cᵢ, ∇◌rᵢ, max_visc = nothing, min_dt_force = nothing; @@ -52,7 +54,7 @@ using LinearAlgebra dρdt_acc = zero(dρdtI[i]) acc_acc = zero(Acceleration[i]) CellIndex = Cells[i] - CellListIndex = get(CellDict, CellIndex, 1) + CellListIndex = CellLookupIndex(CellLookup, CellIndex, 1) SameCellStart = ParticleRanges[CellListIndex] SameCellEnd = ParticleRanges[CellListIndex + 1] - 1 NeighborCellIndices = NeighborCellLists[CellListIndex] @@ -91,7 +93,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 +113,7 @@ using LinearAlgebra kernel_acc = zero(Kernel[i]) kernel_grad_acc = zero(KernelGradient[i]) CellIndex = Cells[i] - CellListIndex = get(CellDict, CellIndex, 1) + CellListIndex = CellLookupIndex(CellLookup, CellIndex, 1) SameCellStart = ParticleRanges[CellListIndex] SameCellEnd = ParticleRanges[CellListIndex + 1] - 1 NeighborCellIndices = NeighborCellLists[CellListIndex] @@ -156,7 +158,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 +177,7 @@ using LinearAlgebra shift_c_acc = zero(∇Cᵢ[i]) shift_r_acc = zero(∇◌rᵢ[i]) CellIndex = Cells[i] - CellListIndex = get(CellDict, CellIndex, 1) + CellListIndex = CellLookupIndex(CellLookup, CellIndex, 1) SameCellStart = ParticleRanges[CellListIndex] SameCellEnd = ParticleRanges[CellListIndex + 1] - 1 NeighborCellIndices = NeighborCellLists[CellListIndex] @@ -220,7 +222,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 +245,7 @@ using LinearAlgebra shift_c_acc = zero(∇Cᵢ[i]) shift_r_acc = zero(∇◌rᵢ[i]) CellIndex = Cells[i] - CellListIndex = get(CellDict, CellIndex, 1) + CellListIndex = CellLookupIndex(CellLookup, CellIndex, 1) SameCellStart = ParticleRanges[CellListIndex] SameCellEnd = ParticleRanges[CellListIndex + 1] - 1 NeighborCellIndices = NeighborCellLists[CellListIndex] @@ -291,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, CellDict, Position, + SimConstants, ParticleRanges, CellLookup, Position, Density, GhostPoints, GhostNormals, ParticleType, bᵧ, Aᵧ; ParticleOrder) where {Dimensions, FloatType, SMode, KMode, BMode, LMode} @@ -314,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 = get(CellDict, SCellIndex, 1) + NeighborIdx = CellLookupIndex(CellLookup, SCellIndex, 1) StartIndex_ = ParticleRanges[NeighborIdx] EndIndex_ = ParticleRanges[NeighborIdx + 1] - 1 @@ -632,7 +634,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 +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, 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 +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, CellDict, + ParticleRanges, UniqueCells, CellLookup, ParticleOrder, CellOffsets, NeighborCellLists, dρdtI, Velocityₙ⁺, Positionₙ⁺, ρₙ⁺, ∇Cᵢ, ∇◌rᵢ, @@ -747,10 +749,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 +762,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 +793,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 +804,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 +869,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 = InitializeCellIndexLookup(Val(Dimensions)) FullStencil = ConstructStencil(Val(Dimensions)) NeighborCellLists = [Int[] for _ in 1:length(UniqueCells)] ParticleOrder = zeros(Int, NumberOfPoints) @@ -909,7 +911,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/SPHExample.jl b/src/SPHExample.jl index a56895ea..a31df067 100644 --- a/src/SPHExample.jl +++ b/src/SPHExample.jl @@ -58,7 +58,8 @@ module SPHExample using .SPHNeighborList export ConstructStencil, ExtractCells!, UpdateNeighbors!, - BuildNeighborCellLists!, ComputeCellParticleCounts, ComputeCellNeighborCounts + BuildNeighborCellLists!, ComputeCellParticleCounts, ComputeCellNeighborCounts, + CellIndexLookup, InitializeCellIndexLookup using .SPHCellList export NeighborLoop!, ComputeInteractions!, RunSimulation diff --git a/src/SPHNeighborList.jl b/src/SPHNeighborList.jl index b94fc07e..5e8d633b 100644 --- a/src/SPHNeighborList.jl +++ b/src/SPHNeighborList.jl @@ -1,15 +1,63 @@ module SPHNeighborList export ConstructStencil, ExtractCells!, UpdateNeighbors!, - BuildNeighborCellLists!, ComputeCellParticleCounts, ComputeCellNeighborCounts, UpdateΔx! + BuildNeighborCellLists!, ComputeCellParticleCounts, ComputeCellNeighborCounts, UpdateΔx!, + CellIndexLookup, InitializeCellIndexLookup, CellLookupIndex using StaticArrays +mutable struct CellIndexLookup{D} + Grid::Vector{Int} + MinCell::MVector{D, Int} + MaxCell::MVector{D, Int} + Strides::MVector{D, Int} +end + +function InitializeCellIndexLookup(::Val{D}) where D + zeros_tuple = ntuple(_ -> 0, D) + return CellIndexLookup{D}( + Int[], + MVector{D, Int}(zeros_tuple), + MVector{D, Int}(zeros_tuple), + MVector{D, Int}(zeros_tuple), + ) +end + +@inline function PrepareCellLookup!(CellLookup::CellIndexLookup{D}) where D + total = 1 + @inbounds for d in 1:D + dim = CellLookup.MaxCell[d] - CellLookup.MinCell[d] + 1 + CellLookup.Strides[d] = total + total *= dim + end + resize!(CellLookup.Grid, total) + fill!(CellLookup.Grid, 0) + return nothing +end + +@inline function CellLinearIndex(CellLookup::CellIndexLookup{D}, Cell::CartesianIndex{D}) where D + idx = 1 + @inbounds for d in 1:D + idx += (Cell[d] - CellLookup.MinCell[d]) * CellLookup.Strides[d] + end + return idx +end + +@inline function CellLookupIndex(CellLookup::CellIndexLookup{D}, Cell::CartesianIndex{D}, Default) where D + @inbounds for d in 1:D + value = Cell[d] + if value < CellLookup.MinCell[d] || value > CellLookup.MaxCell[d] + return Default + end + end + return CellLookup.Grid[CellLinearIndex(CellLookup, Cell)] +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 +75,7 @@ function BuildNeighborCellLists!(NeighborCellLists, FullStencil, UniqueCellsView Cell = UniqueCellsView[CellIndex] for Offset in FullStencil NeighborCell = Cell + Offset - NeighborIndex = get(CellDict, NeighborCell, 0) + NeighborIndex = CellLookupIndex(CellLookup, NeighborCell, 0) if NeighborIndex != 0 && NeighborIndex != CellIndex StartIndex = ParticleRanges[NeighborIndex] EndIndex = ParticleRanges[NeighborIndex + 1] - 1 @@ -66,6 +114,38 @@ end return nothing end +@inline function ExtractCellsAndBounds!(Particles, InverseCutOff, MinCell, MaxCell) + Cells = Particles.Cells + Positions = Particles.Position + if isempty(Cells) + return false + end + + FirstCell = CartesianIndex(map(X -> MapFloor(X, InverseCutOff), Tuple(Positions[1]))) + Cells[1] = FirstCell + @inbounds for d in eachindex(MinCell) + MinCell[d] = FirstCell[d] + MaxCell[d] = FirstCell[d] + end + + if length(Cells) > 1 + @inbounds @simd ivdep for Index in 2:length(Cells) + Cell = CartesianIndex(map(X -> MapFloor(X, InverseCutOff), Tuple(Positions[Index]))) + Cells[Index] = Cell + @inbounds for d in eachindex(MinCell) + value = Cell[d] + if value < MinCell[d] + MinCell[d] = value + elseif value > MaxCell[d] + MaxCell[d] = value + end + end + end + end + + return true +end + """ Updates the neighbor list without sorting particle storage. @@ -73,22 +153,26 @@ 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) - ExtractCells!(Particles, InverseCutOff) + UniqueCells, CellLookup, ParticleOrder, CellOffsets) + HasParticles = ExtractCellsAndBounds!(Particles, InverseCutOff, CellLookup.MinCell, CellLookup.MaxCell) + if !HasParticles + return 0 + end + PrepareCellLookup!(CellLookup) Cells = @views Particles.Cells ParticleRanges[1] = 1 IndexCounter = 1 - empty!(CellDict) fill!(CellOffsets, zero(eltype(CellOffsets))) @inbounds for Index in eachindex(Cells) Cell = Cells[Index] - CellIndex = get(CellDict, Cell, 0) + LinearIndex = CellLinearIndex(CellLookup, Cell) + CellIndex = CellLookup.Grid[LinearIndex] if CellIndex == 0 IndexCounter += 1 CellIndex = IndexCounter - CellDict[Cell] = CellIndex + CellLookup.Grid[LinearIndex] = CellIndex UniqueCells[CellIndex] = Cell end CellOffsets[CellIndex] += 1 @@ -104,7 +188,7 @@ function UpdateNeighbors!(Particles, InverseCutOff, ParticleRanges, ParticleRanges[IndexCounter + 1] = RunningIndex @inbounds for Index in eachindex(Cells) - CellIndex = CellDict[Cells[Index]] + CellIndex = CellLookup.Grid[CellLinearIndex(CellLookup, Cells[Index])] TargetIndex = CellOffsets[CellIndex] ParticleOrder[TargetIndex] = Index CellOffsets[CellIndex] = TargetIndex + 1