diff --git a/src/SPHCellList.jl b/src/SPHCellList.jl index ffd176f4..290e1715 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, DefaultNeighborLipFactor, ExtractCells!, MapFloor, NeighborListInverseCutOff, UpdateNeighbors!, UpdateΔx! using StaticArrays using StructArrays: StructArray, foreachfield @@ -297,7 +297,7 @@ using LinearAlgebra return nothing end - f(SimKernel, GhostPoint) = CartesianIndex(map(x -> MapFloor(x, SimKernel.H⁻¹), Tuple(GhostPoint))) + f(SimKernel, GhostPoint; LipFactor = DefaultNeighborLipFactor) = CartesianIndex(map(x -> MapFloor(x, NeighborListInverseCutOff(SimKernel.H⁻¹, typeof(SimKernel.H⁻¹)(LipFactor))), Tuple(GhostPoint))) function NeighborLoopMDBC!(SimKernel, SimMetaData::SimulationMetaData{Dimensions, FloatType, SMode, KMode, BMode, LMode}, SimConstants, ParticleRanges, CellDict, Position, @@ -724,7 +724,8 @@ using LinearAlgebra @timeit SimMetaData.HourGlass "01 Calculate IndexCounter" begin SimMetaData.Δx = UpdateΔx!(SimMetaData.Δx, Positionₙ⁺, SimParticles.Position) - ShouldRebuild = SimMetaData.Δx >= SimKernel.h + NeighborListLipFactor = FloatType(DefaultNeighborLipFactor) + ShouldRebuild = SimMetaData.Δx >= SimKernel.h * NeighborListLipFactor # println("Δx: ", Δx, "h: ", SimKernel.h," dt: ", SimMetaData.CurrentTimeStep, " Iteration: ", SimMetaData.Iteration, " TotalTime: ", SimMetaData.TotalTime, " OutputIterationCounter: ", SimMetaData.OutputIterationCounter) @@ -735,7 +736,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 - @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⁻¹, SortingScratchSpace, ParticleRanges, UniqueCells, CellDict; LipFactor = NeighborListLipFactor) SimMetaData.Δx = zero(eltype(dρdtI)) UniqueCellsView = view(UniqueCells, 1:SimMetaData.IndexCounter) BuildNeighborCellLists!(NeighborCellLists, FullStencil, UniqueCellsView, ParticleRanges, CellDict) diff --git a/src/SPHNeighborList.jl b/src/SPHNeighborList.jl index 48d245ea..5c804a30 100644 --- a/src/SPHNeighborList.jl +++ b/src/SPHNeighborList.jl @@ -1,9 +1,16 @@ module SPHNeighborList -export ConstructStencil, ExtractCells!, UpdateNeighbors!, BuildNeighborCellLists!, ComputeCellParticleCounts, ComputeCellNeighborCounts, UpdateΔx! +export ConstructStencil, ExtractCells!, UpdateNeighbors!, BuildNeighborCellLists!, ComputeCellParticleCounts, ComputeCellNeighborCounts, UpdateΔx!, NeighborListInverseCutOff using StaticArrays +const DefaultNeighborLipFactor = 1.2 + +@inline function NeighborListInverseCutOff(InverseCutOff::T, LipFactor::T = T(DefaultNeighborLipFactor)) where {T<:Real} + @assert LipFactor > zero(T) "Neighbor list lip factor must be positive" + return InverseCutOff / LipFactor +end + function ConstructStencil(V::Val{d}) where d return CartesianIndices(ntuple(_ -> -1:1, V)) end @@ -42,7 +49,8 @@ Extracts the cells for each particle based on their positions and the inverse cu # Arguments - `Particles`: The particles whose cells are to be extracted. -- `::Val{InverseCutOff}`: The inverse cutoff value used for cell extraction. +- `InverseCutOff`: The inverse cutoff value used for cell extraction. +- `LipFactor`: Scaling factor to expand cell sizes around their centers. # Returns - `nothing`: This function modifies the `Particles` in place. @@ -55,9 +63,10 @@ Extracts the cells for each particle based on their positions and the inverse cu Int(sign(X)) * unsafe_trunc(Int, muladd(abs(X), InverseCutOff, 0.5)) end -@inline function ExtractCells!(Particles, InverseCutOff) +@inline function ExtractCells!(Particles, InverseCutOff, LipFactor = DefaultNeighborLipFactor) + EffectiveInverseCutOff = NeighborListInverseCutOff(InverseCutOff, typeof(InverseCutOff)(LipFactor)) @inbounds @simd ivdep for Index ∈ eachindex(Particles.Cells) - Particles.Cells[Index] = CartesianIndex(map(X -> MapFloor(X, InverseCutOff), Tuple(Particles.Position[Index]))) + Particles.Cells[Index] = CartesianIndex(map(X -> MapFloor(X, EffectiveInverseCutOff), Tuple(Particles.Position[Index]))) end return nothing end @@ -71,13 +80,15 @@ Updates the neighbor list and sorts particles by their cell indices. - `SortingScratchSpace`: Scratch space for sorting. - `ParticleRanges`: Array to store the ranges of particles in each cell. - `UniqueCells`: Array to store the unique cells. +- `LipFactor`: Scaling factor to expand cell sizes around their centers. # Returns - `IndexCounter`: The number of unique cells identified. """ function UpdateNeighbors!(Particles, InverseCutOff, SortingScratchSpace, - ParticleRanges, UniqueCells, CellDict) - ExtractCells!(Particles, InverseCutOff) + ParticleRanges, UniqueCells, CellDict; + LipFactor = DefaultNeighborLipFactor) + ExtractCells!(Particles, InverseCutOff, LipFactor) sort!(Particles, by = p -> p.Cells; scratch=SortingScratchSpace) Cells = @views Particles.Cells