Skip to content

Use cell-index map for neighbor lookups and refactor MDBC neighbor loop #198

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
AhmedSalih3d wants to merge 3 commits into
base: BeforeLidDrivenCavityExample
Choose a base branch
from
+106 −45
Open

Use cell-index map for neighbor lookups and refactor MDBC neighbor loop #198

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
f375d87
Refactor MDBC neighbor lookup
AhmedSalih3d Jul 12, 2026
43d114c
Flatten neighbor cell lists
AhmedSalih3d Jul 12, 2026
d150722
Merge branch 'codex/optimize-neighborloop-algorithm' of https://githu…
AhmedSalih3d Jul 12, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
95 changes: 51 additions & 44 deletions src/SPHCellList.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -15,7 +15,7 @@ using ..OpenExternalPrograms
using ..SPHKernels
using ..SPHViscosityModels
using ..SPHDensityDiffusionModels
using ..SPHNeighborList: BuildNeighborCellLists!, ComputeCellNeighborCounts, ComputeCellParticleCounts, ConstructStencil, ExtractCells!, FindCellIndex, MapFloor, UpdateNeighbors!, UpdateΔx!
using ..SPHNeighborList: BuildNeighborCellLists!, ComputeCellNeighborCounts, ComputeCellParticleCounts, ConstructStencil, ExtractCells!, MapFloor, UpdateNeighbors!, UpdateΔx!

using Base.Threads: @threads
using Bumper: @alloc, @no_escape
Expand All @@ -29,7 +29,7 @@ using TimerOutputs: @timeit
function NeighborLoopPerParticle!(SimDensityDiffusion::SDD, SimViscosity::SV, SimKernel,
SimMetaData::SimulationMetaData{D,T,NoShifting,NoKernelOutput,B,L},
SimConstants, SimParticles, ParticleRanges,
CellListIndices, NeighborCellLists, dρdtI,
CellListIndices, NeighborCellListOffsets, NeighborCellListIndices, dρdtI,
Acceleration, ∇Cᵢ,
∇◌rᵢ, AccelerationMax;
Position = SimParticles.Position,
Expand All @@ -46,7 +46,8 @@ using TimerOutputs: @timeit
CellListIndex = CellListIndices[i]
SameCellStart = ParticleRanges[CellListIndex]
SameCellEnd = ParticleRanges[CellListIndex + 1] - 1
NeighborCellIndices = NeighborCellLists[CellListIndex]
NeighborCellStart = NeighborCellListOffsets[CellListIndex]
NeighborCellEnd = NeighborCellListOffsets[CellListIndex + 1] - 1

@inbounds for j in SameCellStart:(i - 1)
dρdt_acc, acc_acc = ComputeInteractionsPerParticle!(
Expand All @@ -62,7 +63,8 @@ using TimerOutputs: @timeit
Velocity, ParticleType, dρdt_acc, acc_acc, i, j,
)
end
for NeighborIdx in NeighborCellIndices
for NeighborCellCursor in NeighborCellStart:NeighborCellEnd
NeighborIdx = NeighborCellListIndices[NeighborCellCursor]
StartIndex_ = ParticleRanges[NeighborIdx]
EndIndex_ = ParticleRanges[NeighborIdx + 1] - 1
@inbounds for j in StartIndex_:EndIndex_
Expand All @@ -85,7 +87,7 @@ using TimerOutputs: @timeit
function NeighborLoopPerParticle!(SimDensityDiffusion::SDD, SimViscosity::SV, SimKernel,
SimMetaData::SimulationMetaData{D,T,NoShifting,K,B,L},
SimConstants, SimParticles, ParticleRanges,
CellListIndices, NeighborCellLists, dρdtI,
CellListIndices, NeighborCellListOffsets, NeighborCellListIndices, dρdtI,
Acceleration, ∇Cᵢ,
∇◌rᵢ, AccelerationMax;
Position = SimParticles.Position,
Expand All @@ -106,7 +108,8 @@ using TimerOutputs: @timeit
CellListIndex = CellListIndices[i]
SameCellStart = ParticleRanges[CellListIndex]
SameCellEnd = ParticleRanges[CellListIndex + 1] - 1
NeighborCellIndices = NeighborCellLists[CellListIndex]
NeighborCellStart = NeighborCellListOffsets[CellListIndex]
NeighborCellEnd = NeighborCellListOffsets[CellListIndex + 1] - 1

@inbounds for j in SameCellStart:(i - 1)
dρdt_acc, acc_acc, kernel_acc, kernel_grad_acc =
Expand All @@ -126,7 +129,8 @@ using TimerOutputs: @timeit
kernel_grad_acc, i, j,
)
end
for NeighborIdx in NeighborCellIndices
for NeighborCellCursor in NeighborCellStart:NeighborCellEnd
NeighborIdx = NeighborCellListIndices[NeighborCellCursor]
StartIndex_ = ParticleRanges[NeighborIdx]
EndIndex_ = ParticleRanges[NeighborIdx + 1] - 1
@inbounds for j in StartIndex_:EndIndex_
Expand All @@ -153,7 +157,7 @@ using TimerOutputs: @timeit
function NeighborLoopPerParticle!(SimDensityDiffusion::SDD, SimViscosity::SV, SimKernel,
SimMetaData::SimulationMetaData{D,T,S,NoKernelOutput,B,L},
SimConstants, SimParticles, ParticleRanges,
CellListIndices, NeighborCellLists, dρdtI,
CellListIndices, NeighborCellListOffsets, NeighborCellListIndices, dρdtI,
Acceleration, ∇Cᵢ,
∇◌rᵢ, AccelerationMax;
Position = SimParticles.Position,
Expand All @@ -172,7 +176,8 @@ using TimerOutputs: @timeit
CellListIndex = CellListIndices[i]
SameCellStart = ParticleRanges[CellListIndex]
SameCellEnd = ParticleRanges[CellListIndex + 1] - 1
NeighborCellIndices = NeighborCellLists[CellListIndex]
NeighborCellStart = NeighborCellListOffsets[CellListIndex]
NeighborCellEnd = NeighborCellListOffsets[CellListIndex + 1] - 1

@inbounds for j in SameCellStart:(i - 1)
dρdt_acc, acc_acc, shift_c_acc, shift_r_acc =
Expand All @@ -192,7 +197,8 @@ using TimerOutputs: @timeit
shift_r_acc, i, j,
)
end
for NeighborIdx in NeighborCellIndices
for NeighborCellCursor in NeighborCellStart:NeighborCellEnd
NeighborIdx = NeighborCellListIndices[NeighborCellCursor]
StartIndex_ = ParticleRanges[NeighborIdx]
EndIndex_ = ParticleRanges[NeighborIdx + 1] - 1
@inbounds for j in StartIndex_:EndIndex_
Expand All @@ -219,7 +225,7 @@ using TimerOutputs: @timeit
function NeighborLoopPerParticle!(SimDensityDiffusion::SDD, SimViscosity::SV, SimKernel,
SimMetaData::SimulationMetaData{D,T,S,K,B,L},
SimConstants, SimParticles, ParticleRanges,
CellListIndices, NeighborCellLists, dρdtI,
CellListIndices, NeighborCellListOffsets, NeighborCellListIndices, dρdtI,
Acceleration, ∇Cᵢ,
∇◌rᵢ, AccelerationMax;
Position = SimParticles.Position,
Expand All @@ -243,7 +249,8 @@ using TimerOutputs: @timeit
CellListIndex = CellListIndices[i]
SameCellStart = ParticleRanges[CellListIndex]
SameCellEnd = ParticleRanges[CellListIndex + 1] - 1
NeighborCellIndices = NeighborCellLists[CellListIndex]
NeighborCellStart = NeighborCellListOffsets[CellListIndex]
NeighborCellEnd = NeighborCellListOffsets[CellListIndex + 1] - 1

@inbounds for j in SameCellStart:(i - 1)
dρdt_acc, acc_acc, kernel_acc, kernel_grad_acc, shift_c_acc,
Expand All @@ -263,7 +270,8 @@ using TimerOutputs: @timeit
kernel_grad_acc, shift_c_acc, shift_r_acc, i, j,
)
end
for NeighborIdx in NeighborCellIndices
for NeighborCellCursor in NeighborCellStart:NeighborCellEnd
NeighborIdx = NeighborCellListIndices[NeighborCellCursor]
StartIndex_ = ParticleRanges[NeighborIdx]
EndIndex_ = ParticleRanges[NeighborIdx + 1] - 1
@inbounds for j in StartIndex_:EndIndex_
Expand All @@ -289,17 +297,15 @@ using TimerOutputs: @timeit
return nothing
end

f(SimKernel, GhostPoint) = CartesianIndex(map(x -> MapFloor(x, SimKernel.H⁻¹), Tuple(GhostPoint)))
@inline GhostCellIndex(SimKernel, GhostPoint) = CartesianIndex(map(x -> MapFloor(x, SimKernel.H⁻¹), Tuple(GhostPoint)))
function NeighborLoopMDBC!(SimKernel,
SimMetaData::SimulationMetaData{Dimensions, FloatType, SMode, KMode, BMode, LMode},
SimConstants, ParticleRanges, UniqueCellsView,
SimParticles, bᵧ, Aᵧ) where {Dimensions, FloatType, SMode, KMode, BMode, LMode}
SimConstants, FullStencil, ParticleRanges,
CellIndexMap, SimParticles, bᵧ, Aᵧ) where {Dimensions, FloatType, SMode, KMode, BMode, LMode}

@unpack Position, Density, GhostPoints, GhostNormals = SimParticles
@unpack Position, Density, GhostPoints = SimParticles
ParticleType = SimParticles.Type

FullStencil = ConstructStencil(Val(Dimensions))

@inbounds @threads for iter in eachindex(GhostPoints)
GhostPoint = GhostPoints[iter]

Expand All @@ -309,14 +315,12 @@ using TimerOutputs: @timeit
A_acc = zero(Aᵧ[iter]) # an SMatrix{D+1,D+1,FloatType}

# compute and accumulate into the locals
GhostCellIndex = f(SimKernel, GhostPoints[iter])
@inbounds for offset ∈ FullStencil
SCellIndex = GhostCellIndex + offset

# Returns a range, x>:x for exact match and x=:x for no match
# utilizes that it is a sorted array and requires no isequal constructor,
# so I prefer this for now
NeighborIdx = FindCellIndex(UniqueCellsView, SCellIndex)
CellIndex = GhostCellIndex(SimKernel, GhostPoint)
@inbounds for Offset ∈ FullStencil
NeighborIdx = get(CellIndexMap, CellIndex + Offset, 0)
if iszero(NeighborIdx)
continue
end

StartIndex_ = ParticleRanges[NeighborIdx]
EndIndex_ = ParticleRanges[NeighborIdx + 1] - 1
Expand Down Expand Up @@ -628,14 +632,13 @@ using TimerOutputs: @timeit

function ApplyMDBCBeforeHalf!(SimMetaData::SimulationMetaData{D,T,S,K,SimpleMDBC,L},
SimKernel, SimConstants, SimParticles,
ParticleRanges, UniqueCells
FullStencil, ParticleRanges, CellIndexMap
) where {D,T,S<:ShiftingMode,K<:KernelOutputMode,L<:LogMode}
@no_escape begin
DimensionsPlus = D + 1
bᵧ = @alloc(SVector{DimensionsPlus, T}, length(SimParticles.Position))
Aᵧ = @alloc(SMatrix{DimensionsPlus, DimensionsPlus, T, DimensionsPlus*DimensionsPlus}, length(SimParticles.Position))
UniqueCellsView = view(UniqueCells, 1:SimMetaData.IndexCounter)
NeighborLoopMDBC!(SimKernel, SimMetaData, SimConstants, ParticleRanges, UniqueCellsView, SimParticles, bᵧ, Aᵧ)
NeighborLoopMDBC!(SimKernel, SimMetaData, SimConstants, FullStencil, ParticleRanges, CellIndexMap, SimParticles, bᵧ, Aᵧ)
ApplyMDBCCorrection(SimConstants, SimParticles, bᵧ, Aᵧ)
end

Expand Down Expand Up @@ -696,7 +699,7 @@ using TimerOutputs: @timeit
SimConstants, SimParticles, FullStencil,
ParticleRanges, UniqueCells, CellListIndices,
SortingScratchSpace,
NeighborCellLists, dρdtI, Velocityₙ⁺,
NeighborCellListOffsets, NeighborCellListIndices, NeighborCellIndexMap, dρdtI, Velocityₙ⁺,
Positionₙ⁺, ρₙ⁺, ∇Cᵢ, ∇◌rᵢ,
MotionDefinition::Union{
Nothing,
Expand Down Expand Up @@ -728,15 +731,15 @@ using TimerOutputs: @timeit

SimMetaData.IndexCounter = UpdateNeighbors!(SimParticles, SimKernel.H⁻¹, SortingScratchSpace, ParticleRanges, UniqueCells, CellListIndices)
UniqueCellsView = view(UniqueCells, 1:SimMetaData.IndexCounter)
BuildNeighborCellLists!(NeighborCellLists, FullStencil, UniqueCellsView, ParticleRanges)
BuildNeighborCellLists!(NeighborCellListOffsets, NeighborCellListIndices, FullStencil, UniqueCellsView, ParticleRanges, NeighborCellIndexMap)

if TimeSteppingMode isa SingleNeighborTimeStepping
@timeit SimMetaData.HourGlass "00 Init Pressure" Pressure!(SimParticles.Pressure, SimParticles.Density, SimConstants)
@timeit SimMetaData.HourGlass "00a Init MDBC" ApplyMDBCBeforeHalf!(SimMetaData, SimKernel, SimConstants, SimParticles, ParticleRanges, UniqueCells)
@timeit SimMetaData.HourGlass "00a Init MDBC" ApplyMDBCBeforeHalf!(SimMetaData, SimKernel, SimConstants, SimParticles, FullStencil, ParticleRanges, NeighborCellIndexMap)
@timeit SimMetaData.HourGlass "00b Init NeighborLoop" NeighborLoopPerParticle!(
SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData,
SimConstants, SimParticles, ParticleRanges, CellListIndices,
NeighborCellLists, dρdtI, SimParticles.Acceleration, ∇Cᵢ, ∇◌rᵢ, AccelerationMax,
NeighborCellListOffsets, NeighborCellListIndices, dρdtI, SimParticles.Acceleration, ∇Cᵢ, ∇◌rᵢ, AccelerationMax,
)
end

Expand All @@ -759,20 +762,20 @@ using TimerOutputs: @timeit
@timeit SimMetaData.HourGlass "01a Actual Calculate IndexCounter" SimMetaData.IndexCounter = UpdateNeighbors!(SimParticles, SimKernel.H⁻¹, SortingScratchSpace, ParticleRanges, UniqueCells, CellListIndices)
SimMetaData.Δx = zero(eltype(dρdtI))
UniqueCellsView = view(UniqueCells, 1:SimMetaData.IndexCounter)
BuildNeighborCellLists!(NeighborCellLists, FullStencil, UniqueCellsView, ParticleRanges)
BuildNeighborCellLists!(NeighborCellListOffsets, NeighborCellListIndices, FullStencil, UniqueCellsView, ParticleRanges, NeighborCellIndexMap)
end
end

@timeit SimMetaData.HourGlass "Motion" ProgressMotion(SimParticles, dt₂, MotionDefinition, SimMetaData)

if TimeSteppingMode isa SymplecticTimeStepping
@timeit SimMetaData.HourGlass "02 Pressure" Pressure!(SimParticles.Pressure, SimParticles.Density, SimConstants)
@timeit SimMetaData.HourGlass "03 Apply MDBC before Half TimeStep" ApplyMDBCBeforeHalf!(SimMetaData, SimKernel, SimConstants, SimParticles, ParticleRanges, UniqueCells)
@timeit SimMetaData.HourGlass "03 Apply MDBC before Half TimeStep" ApplyMDBCBeforeHalf!(SimMetaData, SimKernel, SimConstants, SimParticles, FullStencil, ParticleRanges, NeighborCellIndexMap)

@timeit SimMetaData.HourGlass "04 First NeighborLoop" NeighborLoopPerParticle!(
SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData,
SimConstants, SimParticles, ParticleRanges, CellListIndices,
NeighborCellLists, dρdtI, SimParticles.Acceleration, ∇Cᵢ, ∇◌rᵢ, AccelerationMax,
NeighborCellListOffsets, NeighborCellListIndices, dρdtI, SimParticles.Acceleration, ∇Cᵢ, ∇◌rᵢ, AccelerationMax,
)

@timeit SimMetaData.HourGlass "05 Update To Half TimeStep" HalfTimeStep(SimMetaData, SimConstants, SimParticles, Positionₙ⁺, Velocityₙ⁺, ρₙ⁺, dρdtI, dt₂)
Expand All @@ -785,13 +788,13 @@ using TimerOutputs: @timeit
@timeit SimMetaData.HourGlass "08 Second NeighborLoop" NeighborLoopPerParticle!(
SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData,
SimConstants, SimParticles, ParticleRanges, CellListIndices,
NeighborCellLists, dρdtI, SimParticles.Acceleration, ∇Cᵢ, ∇◌rᵢ, AccelerationMax,
NeighborCellListOffsets, NeighborCellListIndices, dρdtI, SimParticles.Acceleration, ∇Cᵢ, ∇◌rᵢ, AccelerationMax,
Position = Positionₙ⁺,
Density = ρₙ⁺,
Velocity = Velocityₙ⁺,
)
else
@timeit SimMetaData.HourGlass "02 Apply MDBC before Half TimeStep" ApplyMDBCBeforeHalf!(SimMetaData, SimKernel, SimConstants, SimParticles, ParticleRanges, UniqueCells)
@timeit SimMetaData.HourGlass "02 Apply MDBC before Half TimeStep" ApplyMDBCBeforeHalf!(SimMetaData, SimKernel, SimConstants, SimParticles, FullStencil, ParticleRanges, NeighborCellIndexMap)

@timeit SimMetaData.HourGlass "03 Update To Half TimeStep" HalfTimeStep(SimMetaData, SimConstants, SimParticles, Positionₙ⁺, Velocityₙ⁺, ρₙ⁺, dρdtI, dt₂)

Expand All @@ -803,7 +806,7 @@ using TimerOutputs: @timeit
@timeit SimMetaData.HourGlass "06 NeighborLoop" NeighborLoopPerParticle!(
SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData,
SimConstants, SimParticles, ParticleRanges, CellListIndices,
NeighborCellLists, dρdtI, SimParticles.Acceleration, ∇Cᵢ, ∇◌rᵢ, AccelerationMax,
NeighborCellListOffsets, NeighborCellListIndices, dρdtI, SimParticles.Acceleration, ∇Cᵢ, ∇◌rᵢ, AccelerationMax,
Position = Positionₙ⁺,
Density = ρₙ⁺,
Velocity = Velocityₙ⁺,
Expand Down Expand Up @@ -911,7 +914,9 @@ using TimerOutputs: @timeit
UniqueCells = zeros(CartesianIndex{Dimensions}, NumberOfPoints)
CellListIndices = zeros(Int, NumberOfPoints)
FullStencil = ConstructStencil(Val(Dimensions))
NeighborCellLists = [Int[] for _ in 1:length(UniqueCells)]
NeighborCellListOffsets = zeros(Int, length(UniqueCells) + 1)
NeighborCellListIndices = Int[]
NeighborCellIndexMap = Dict{CartesianIndex{Dimensions}, Int}()
_, SortingScratchSpace = Base.Sort.make_scratch(nothing, eltype(SimParticles), NumberOfPoints)

output = SetupVTKOutput(SimMetaData, SimParticles, SimKernel, Dimensions)
Expand All @@ -930,7 +935,8 @@ using TimerOutputs: @timeit
)
cell_neighbor_counts = ComputeCellNeighborCounts(
ParticleRanges,
NeighborCellLists,
NeighborCellListOffsets,
NeighborCellListIndices,
SimMetaData.IndexCounter,
)
end
Expand All @@ -952,7 +958,7 @@ using TimerOutputs: @timeit
SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData,
SimConstants, SimParticles, FullStencil, ParticleRanges,
UniqueCells, CellListIndices, SortingScratchSpace,
NeighborCellLists, dρdtI, Velocityₙ⁺, Positionₙ⁺, ρₙ⁺,
NeighborCellListOffsets, NeighborCellListIndices, NeighborCellIndexMap, dρdtI, Velocityₙ⁺, Positionₙ⁺, ρₙ⁺,
∇Cᵢ, ∇◌rᵢ, MotionDefinition,
)
push!(SimMetaData.TimeSteps, SimMetaData.CurrentTimeStep)
Expand All @@ -971,7 +977,8 @@ using TimerOutputs: @timeit
)
cell_neighbor_counts = ComputeCellNeighborCounts(
ParticleRanges,
NeighborCellLists,
NeighborCellListOffsets,
NeighborCellListIndices,
length(UniqueCellsView),
)
end
Expand Down
Loading