diff --git a/README.md b/README.md index e1bf6bb4..c23bc265 100644 --- a/README.md +++ b/README.md @@ -88,6 +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. +Neighbor rebuilds use a per-cell ordering buffer without reordering particle +storage. ## Help diff --git a/src/SPHCellList.jl b/src/SPHCellList.jl index fef4e636..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,14 +37,15 @@ 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; Position = SimParticles.Position, Density = SimParticles.Density, Pressure = SimParticles.Pressure, - Velocity = SimParticles.Velocity) where {D,T, + Velocity = SimParticles.Velocity, + ParticleOrder) where {D,T, B<:MDBCMode,L<:LogMode, SDD<:SPHDensityDiffusion, SV<:SPHViscosity} @@ -51,33 +54,30 @@ 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] - @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, - ) + @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, jIndex, + ) + end 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 end @@ -93,14 +93,15 @@ 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; Position = SimParticles.Position, Density = SimParticles.Density, Pressure = SimParticles.Pressure, - Velocity = SimParticles.Velocity) where {D,T, + Velocity = SimParticles.Velocity, + ParticleOrder) where {D,T, K<:KernelOutputMode, B<:MDBCMode,L<:LogMode, SDD<:SPHDensityDiffusion, @@ -112,39 +113,34 @@ 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] - @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, - ) + @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, jIndex, + ) + end 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 end @@ -162,14 +158,15 @@ 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; Position = SimParticles.Position, Density = SimParticles.Density, Pressure = SimParticles.Pressure, - Velocity = SimParticles.Velocity) where {D,T, + Velocity = SimParticles.Velocity, + ParticleOrder) where {D,T, S<:ShiftingMode,B<:MDBCMode, L<:LogMode,SDD<:SPHDensityDiffusion, SV<:SPHViscosity} @@ -180,39 +177,34 @@ 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] - @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, - ) + @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, jIndex, + ) + end 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 end @@ -230,14 +222,15 @@ 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; Position = SimParticles.Position, Density = SimParticles.Density, Pressure = SimParticles.Pressure, - Velocity = SimParticles.Velocity) where {D,T, + Velocity = SimParticles.Velocity, + ParticleOrder) where {D,T, S<:ShiftingMode, K<:KernelOutputMode, B<:MDBCMode,L<:LogMode, @@ -252,39 +245,34 @@ 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] - @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, - ) + @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, jIndex, + ) + end 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 end @@ -305,9 +293,10 @@ 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ᵧ) where {Dimensions, FloatType, SMode, KMode, BMode, LMode} + bᵧ, Aᵧ; + ParticleOrder) where {Dimensions, FloatType, SMode, KMode, BMode, LMode} FullStencil = ConstructStencil(Val(Dimensions)) @@ -327,16 +316,17 @@ 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 - for j in StartIndex_:EndIndex_ + @inbounds for j in StartIndex_:EndIndex_ # change ComputeInteractions to take & return contributions, e.g.: + jIndex = ParticleOrder[j] bΔ, AΔ = ComputeInteractionsMDBC!(SimKernel, SimMetaData, SimConstants, Position, Density, ParticleType, - GhostPoints, iter, j) + GhostPoints, iter, jIndex) b_acc += bΔ A_acc += AΔ end @@ -644,14 +634,15 @@ using LinearAlgebra function ApplyMDBCBeforeHalf!(SimMetaData::SimulationMetaData{D,T,S,K,SimpleMDBC,L}, SimKernel, SimConstants, SimParticles, - ParticleRanges, CellDict, Position, Density, - GhostPoints, GhostNormals, ParticleType + ParticleRanges, CellLookup, Position, Density, + GhostPoints, GhostNormals, ParticleType; + ParticleOrder ) 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, CellLookup, Position, Density, GhostPoints, GhostNormals, ParticleType, bᵧ, Aᵧ; ParticleOrder = ParticleOrder) ApplyMDBCCorrection(SimConstants, SimParticles, bᵧ, Aᵧ) end @@ -708,8 +699,8 @@ using LinearAlgebra @inbounds function SimulationLoop(SimDensityDiffusion::SDD, SimViscosity::SV, SimKernel, SimMetaData::SimulationMetaData{Dimensions, FloatType, SMode, KMode, BMode, LMode}, SimConstants, SimParticles, FullStencil, - ParticleRanges, UniqueCells, CellDict, - SortingScratchSpace, + ParticleRanges, UniqueCells, CellLookup, + ParticleOrder, CellOffsets, NeighborCellLists, dρdtI, Velocityₙ⁺, Positionₙ⁺, ρₙ⁺, ∇Cᵢ, ∇◌rᵢ, MotionDefinition::Union{ @@ -758,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⁻¹, SortingScratchSpace, ParticleRanges, UniqueCells, CellDict) + @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 @@ -771,20 +762,22 @@ 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) + @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, ) end @@ -800,9 +793,10 @@ 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, Position = Positionₙ⁺, Density = ρₙ⁺, Velocity = Velocityₙ⁺, @@ -810,9 +804,10 @@ 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, Position = Positionₙ⁺, Density = ρₙ⁺, Velocity = Velocityₙ⁺, @@ -833,6 +828,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 @@ -859,10 +869,11 @@ 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)] - _, SortingScratchSpace = Base.Sort.make_scratch(nothing, eltype(SimParticles), NumberOfPoints) + ParticleOrder = zeros(Int, NumberOfPoints) + CellOffsets = zeros(Int, length(ParticleRanges)) output = SetupVTKOutput(SimMetaData, SimParticles, SimKernel, Dimensions) @@ -900,7 +911,7 @@ using LinearAlgebra @timeit SimMetaData.HourGlass "00 SimulationLoop" SimulationLoop( SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData, SimConstants, SimParticles, FullStencil, ParticleRanges, - UniqueCells, CellDict, SortingScratchSpace, + UniqueCells, CellLookup, ParticleOrder, CellOffsets, NeighborCellLists, dρdtI, Velocityₙ⁺, Positionₙ⁺, ρₙ⁺, ∇Cᵢ, ∇◌rᵢ, MotionDefinition, ) @@ -911,24 +922,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 @@ -938,11 +934,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 diff --git a/src/SPHExample.jl b/src/SPHExample.jl index 28ef3ae7..a31df067 100644 --- a/src/SPHExample.jl +++ b/src/SPHExample.jl @@ -57,7 +57,9 @@ module SPHExample export SimulationConstants using .SPHNeighborList - export ConstructStencil, ExtractCells!, UpdateNeighbors!, BuildNeighborCellLists!, ComputeCellParticleCounts, ComputeCellNeighborCounts + export ConstructStencil, ExtractCells!, UpdateNeighbors!, + BuildNeighborCellLists!, ComputeCellParticleCounts, ComputeCellNeighborCounts, + CellIndexLookup, InitializeCellIndexLookup using .SPHCellList export NeighborLoop!, ComputeInteractions!, RunSimulation diff --git a/src/SPHNeighborList.jl b/src/SPHNeighborList.jl index 48d245ea..5e8d633b 100644 --- a/src/SPHNeighborList.jl +++ b/src/SPHNeighborList.jl @@ -1,14 +1,63 @@ module SPHNeighborList -export ConstructStencil, ExtractCells!, UpdateNeighbors!, BuildNeighborCellLists!, ComputeCellParticleCounts, ComputeCellNeighborCounts, UpdateΔx! +export ConstructStencil, ExtractCells!, UpdateNeighbors!, + 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) @@ -22,14 +71,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 = CellLookupIndex(CellLookup, 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 @@ -62,42 +114,85 @@ end return nothing end -""" -Updates the neighbor list and sorts particles by their cell indices. +@inline function ExtractCellsAndBounds!(Particles, InverseCutOff, MinCell, MaxCell) + Cells = Particles.Cells + Positions = Particles.Position + if isempty(Cells) + return false + end -# 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. + 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 -# Returns -- `IndexCounter`: The number of unique cells identified. + 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. + +This builds a per-cell particle ordering buffer so cell ranges can be iterated +without reordering the particle arrays. """ -function UpdateNeighbors!(Particles, InverseCutOff, SortingScratchSpace, - ParticleRanges, UniqueCells, CellDict) - ExtractCells!(Particles, InverseCutOff) +function UpdateNeighbors!(Particles, InverseCutOff, ParticleRanges, + UniqueCells, CellLookup, ParticleOrder, CellOffsets) + HasParticles = ExtractCellsAndBounds!(Particles, InverseCutOff, CellLookup.MinCell, CellLookup.MaxCell) + if !HasParticles + return 0 + end + PrepareCellLookup!(CellLookup) - sort!(Particles, by = p -> p.Cells; scratch=SortingScratchSpace) Cells = @views Particles.Cells - @. 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 + IndexCounter = 1 + fill!(CellOffsets, zero(eltype(CellOffsets))) + + @inbounds for Index in eachindex(Cells) + Cell = Cells[Index] + LinearIndex = CellLinearIndex(CellLookup, Cell) + CellIndex = CellLookup.Grid[LinearIndex] + if CellIndex == 0 + IndexCounter += 1 + CellIndex = IndexCounter + CellLookup.Grid[LinearIndex] = CellIndex + UniqueCells[CellIndex] = Cell end + CellOffsets[CellIndex] += 1 + end + + RunningIndex = 1 + @inbounds for CellIndex in 2:IndexCounter + Count = CellOffsets[CellIndex] + ParticleRanges[CellIndex] = RunningIndex + CellOffsets[CellIndex] = RunningIndex + RunningIndex += Count + end + ParticleRanges[IndexCounter + 1] = RunningIndex + + @inbounds for Index in eachindex(Cells) + CellIndex = CellLookup.Grid[CellLinearIndex(CellLookup, Cells[Index])] + TargetIndex = CellOffsets[CellIndex] + ParticleOrder[TargetIndex] = Index + CellOffsets[CellIndex] = TargetIndex + 1 end - ParticleRanges[IndexCounter + 1] = length(ParticleRanges) return IndexCounter end