From d1424b3bd966f7613ac0e01602886eb9a38628ef Mon Sep 17 00:00:00 2001 From: AhmedSalih3d <36305327+AhmedSalih3d@users.noreply.github.com> Date: Sat, 17 Jan 2026 22:03:18 +0100 Subject: [PATCH 1/8] Localize support buffers in loop --- src/SPHCellList.jl | 25 ++++++++++++++++++------- 1 file changed, 18 insertions(+), 7 deletions(-) diff --git a/src/SPHCellList.jl b/src/SPHCellList.jl index eb55a820..e4c65296 100644 --- a/src/SPHCellList.jl +++ b/src/SPHCellList.jl @@ -688,9 +688,7 @@ using LinearAlgebra SimMetaData::SimulationMetaData{Dimensions, FloatType, SMode, KMode, BMode, LMode}, SimConstants, SimParticles, FullStencil, ParticleRanges, UniqueCells, CellDict, - SortingScratchSpace, - NeighborCellLists, dρdtI, Velocityₙ⁺, - Positionₙ⁺, ρₙ⁺, ∇Cᵢ, ∇◌rᵢ, + SortingScratchSpace, NeighborCellLists, MotionDefinition::Union{ Nothing, AbstractVector{ @@ -718,7 +716,23 @@ using LinearAlgebra @no_escape begin AccelerationMax = @alloc(FloatType, length(Position)) + dρdtI = @alloc(FloatType, length(Position)) + PositionType = eltype(Position) + PositionUnderlyingType = eltype(PositionType) + Velocityₙ⁺ = @alloc(PositionType, length(Position)) + Positionₙ⁺ = @alloc(PositionType, length(Position)) + ρₙ⁺ = @alloc(PositionUnderlyingType, length(Position)) + ∇Cᵢ = SMode <: NoShifting ? Vector{PositionType}(undef, 0) : @alloc(PositionType, length(Position)) + ∇◌rᵢ = SMode <: NoShifting ? Vector{PositionUnderlyingType}(undef, 0) : @alloc(PositionUnderlyingType, length(Position)) dt₂ = dt * 0.5 + fill!(dρdtI, zero(FloatType)) + if !(SMode <: NoShifting) + fill!(∇Cᵢ, zero(PositionType)) + fill!(∇◌rᵢ, zero(PositionUnderlyingType)) + end + copyto!(Positionₙ⁺, Position) + copyto!(Velocityₙ⁺, Velocity) + copyto!(ρₙ⁺, Density) while SimMetaData.TotalTime <= next_output_time(SimMetaData) @timeit SimMetaData.HourGlass "01 Calculate IndexCounter" begin @@ -805,8 +819,6 @@ using LinearAlgebra ) where {Dimensions,FloatType,SMode,KMode,BMode,LMode,SV<:SPHViscosity,SDD<:SPHDensityDiffusion} NumberOfPoints = length(SimParticles) - - dρdtI, Velocityₙ⁺, Positionₙ⁺, ρₙ⁺, ∇Cᵢ, ∇◌rᵢ = AllocateSupportDataStructures(SimMetaData, SimParticles.Position) LoadMDBCNormals!(SimMetaData, SimParticles, ParticleNormalsPath) @@ -859,8 +871,7 @@ using LinearAlgebra SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData, SimConstants, SimParticles, FullStencil, ParticleRanges, UniqueCells, CellDict, SortingScratchSpace, - NeighborCellLists, dρdtI, Velocityₙ⁺, Positionₙ⁺, ρₙ⁺, - ∇Cᵢ, ∇◌rᵢ, MotionDefinition, + NeighborCellLists, MotionDefinition, ) push!(SimMetaData.TimeSteps, SimMetaData.CurrentTimeStep) From d381059e5b62c9aa23cbed1580e737f960e4d3fd Mon Sep 17 00:00:00 2001 From: AhmedSalih3d <36305327+AhmedSalih3d@users.noreply.github.com> Date: Sat, 17 Jan 2026 22:08:07 +0100 Subject: [PATCH 2/8] Use standard arrays in loop --- src/SPHCellList.jl | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/src/SPHCellList.jl b/src/SPHCellList.jl index e4c65296..e64bb46e 100644 --- a/src/SPHCellList.jl +++ b/src/SPHCellList.jl @@ -715,15 +715,15 @@ using LinearAlgebra dt = SimConstants.CFL * (SimKernel.h / SimConstants.c₀) @no_escape begin - AccelerationMax = @alloc(FloatType, length(Position)) - dρdtI = @alloc(FloatType, length(Position)) PositionType = eltype(Position) PositionUnderlyingType = eltype(PositionType) - Velocityₙ⁺ = @alloc(PositionType, length(Position)) - Positionₙ⁺ = @alloc(PositionType, length(Position)) - ρₙ⁺ = @alloc(PositionUnderlyingType, length(Position)) - ∇Cᵢ = SMode <: NoShifting ? Vector{PositionType}(undef, 0) : @alloc(PositionType, length(Position)) - ∇◌rᵢ = SMode <: NoShifting ? Vector{PositionUnderlyingType}(undef, 0) : @alloc(PositionUnderlyingType, length(Position)) + AccelerationMax = similar(Density) + dρdtI = similar(Density) + Velocityₙ⁺ = similar(Position) + Positionₙ⁺ = similar(Position) + ρₙ⁺ = similar(Density) + ∇Cᵢ = SMode <: NoShifting ? Vector{PositionType}(undef, 0) : similar(Position) + ∇◌rᵢ = SMode <: NoShifting ? Vector{PositionUnderlyingType}(undef, 0) : similar(Density) dt₂ = dt * 0.5 fill!(dρdtI, zero(FloatType)) if !(SMode <: NoShifting) From d5a704a85253a878c91d939774a368cc0e2c0345 Mon Sep 17 00:00:00 2001 From: AhmedSalih3d <36305327+AhmedSalih3d@users.noreply.github.com> Date: Sat, 17 Jan 2026 22:49:46 +0100 Subject: [PATCH 3/8] Rebuild neighbors on first step --- src/SPHCellList.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/SPHCellList.jl b/src/SPHCellList.jl index e64bb46e..1fe4623a 100644 --- a/src/SPHCellList.jl +++ b/src/SPHCellList.jl @@ -738,7 +738,7 @@ using LinearAlgebra @timeit SimMetaData.HourGlass "01 Calculate IndexCounter" begin SimMetaData.Δx = UpdateΔx!(SimMetaData.Δx, Positionₙ⁺, SimParticles.Position) - ShouldRebuild = SimMetaData.Δx >= SimKernel.h + ShouldRebuild = SimMetaData.Δx >= SimKernel.h || SimMetaData.IndexCounter == 0 # println("Δx: ", Δx, "h: ", SimKernel.h," dt: ", SimMetaData.CurrentTimeStep, " Iteration: ", SimMetaData.Iteration, " TotalTime: ", SimMetaData.TotalTime, " OutputIterationCounter: ", SimMetaData.OutputIterationCounter) From 2ba07a4882fa3838a128a762aafc245fc68bef52 Mon Sep 17 00:00:00 2001 From: AhmedSalih3d <36305327+AhmedSalih3d@users.noreply.github.com> Date: Sat, 17 Jan 2026 22:55:42 +0100 Subject: [PATCH 4/8] Guard particle range sizing --- src/SPHNeighborList.jl | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/SPHNeighborList.jl b/src/SPHNeighborList.jl index 48d245ea..1e84c7fa 100644 --- a/src/SPHNeighborList.jl +++ b/src/SPHNeighborList.jl @@ -77,6 +77,10 @@ Updates the neighbor list and sorts particles by their cell indices. """ function UpdateNeighbors!(Particles, InverseCutOff, SortingScratchSpace, ParticleRanges, UniqueCells, CellDict) + RequiredLen = length(Particles) + 2 + if length(ParticleRanges) < RequiredLen + resize!(ParticleRanges, RequiredLen) + end ExtractCells!(Particles, InverseCutOff) sort!(Particles, by = p -> p.Cells; scratch=SortingScratchSpace) @@ -97,7 +101,7 @@ function UpdateNeighbors!(Particles, InverseCutOff, SortingScratchSpace, CellDict[Cells[Index]] = IndexCounter end end - ParticleRanges[IndexCounter + 1] = length(ParticleRanges) + ParticleRanges[IndexCounter + 1] = length(Particles) + 1 return IndexCounter end From 247bef76608937bf9dcf9f0c164a05a94e295120 Mon Sep 17 00:00:00 2001 From: AhmedSalih3d <36305327+AhmedSalih3d@users.noreply.github.com> Date: Sat, 17 Jan 2026 23:03:30 +0100 Subject: [PATCH 5/8] Compute half-step fields on demand --- src/SPHCellList.jl | 67 ++++++++++++++++++++++++++++---------- src/SimulationEquations.jl | 22 +++++++++++++ 2 files changed, 72 insertions(+), 17 deletions(-) diff --git a/src/SPHCellList.jl b/src/SPHCellList.jl index 1fe4623a..95d6da91 100644 --- a/src/SPHCellList.jl +++ b/src/SPHCellList.jl @@ -32,6 +32,49 @@ using Base.Threads using LinearAlgebra using Bumper + struct HalfStepPosition{P, V, M, T} + Position::P + Velocity::V + MotionLimiter::M + Δt₂::T + end + + Base.IndexStyle(::Type{<:HalfStepPosition}) = IndexLinear() + Base.axes(half::HalfStepPosition) = axes(half.Position) + Base.size(half::HalfStepPosition) = size(half.Position) + Base.getindex(half::HalfStepPosition, i::Int) = half.Position[i] + half.Velocity[i] * half.Δt₂ * half.MotionLimiter[i] + + struct HalfStepVelocity{V, A, M, T} + Velocity::V + Acceleration::A + MotionLimiter::M + Δt₂::T + end + + Base.IndexStyle(::Type{<:HalfStepVelocity}) = IndexLinear() + Base.axes(half::HalfStepVelocity) = axes(half.Velocity) + Base.size(half::HalfStepVelocity) = size(half.Velocity) + Base.getindex(half::HalfStepVelocity, i::Int) = half.Velocity[i] + half.Acceleration[i] * half.Δt₂ * half.MotionLimiter[i] + + struct HalfStepDensity{D, R, M, T} + Density::D + dρdtI::R + MotionLimiter::M + ρ₀::T + Δt₂::T + end + + Base.IndexStyle(::Type{<:HalfStepDensity}) = IndexLinear() + Base.axes(half::HalfStepDensity) = axes(half.Density) + Base.size(half::HalfStepDensity) = size(half.Density) + function Base.getindex(half::HalfStepDensity, i::Int) + density = half.Density[i] + half.dρdtI[i] * half.Δt₂ + if (density < half.ρ₀) * !Bool(half.MotionLimiter[i]) + density = half.ρ₀ + end + return density + end + function NeighborLoopPerParticle!(SimDensityDiffusion::SDD, SimViscosity::SV, SimKernel, SimMetaData::SimulationMetaData{D,T,NoShifting,NoKernelOutput,B,L}, SimConstants, SimParticles, ParticleRanges, @@ -719,9 +762,6 @@ using LinearAlgebra PositionUnderlyingType = eltype(PositionType) AccelerationMax = similar(Density) dρdtI = similar(Density) - Velocityₙ⁺ = similar(Position) - Positionₙ⁺ = similar(Position) - ρₙ⁺ = similar(Density) ∇Cᵢ = SMode <: NoShifting ? Vector{PositionType}(undef, 0) : similar(Position) ∇◌rᵢ = SMode <: NoShifting ? Vector{PositionUnderlyingType}(undef, 0) : similar(Density) dt₂ = dt * 0.5 @@ -730,14 +770,12 @@ using LinearAlgebra fill!(∇Cᵢ, zero(PositionType)) fill!(∇◌rᵢ, zero(PositionUnderlyingType)) end - copyto!(Positionₙ⁺, Position) - copyto!(Velocityₙ⁺, Velocity) - copyto!(ρₙ⁺, Density) while SimMetaData.TotalTime <= next_output_time(SimMetaData) @timeit SimMetaData.HourGlass "01 Calculate IndexCounter" begin - SimMetaData.Δx = UpdateΔx!(SimMetaData.Δx, Positionₙ⁺, SimParticles.Position) + PositionHalfStep = HalfStepPosition(Position, Velocity, MotionLimiter, dt₂) + SimMetaData.Δx = UpdateΔx!(SimMetaData.Δx, PositionHalfStep, SimParticles.Position) ShouldRebuild = SimMetaData.Δx >= SimKernel.h || SimMetaData.IndexCounter == 0 # println("Δx: ", Δx, "h: ", SimKernel.h," dt: ", SimMetaData.CurrentTimeStep, " Iteration: ", SimMetaData.Iteration, " TotalTime: ", SimMetaData.TotalTime, " OutputIterationCounter: ", SimMetaData.OutputIterationCounter) @@ -771,26 +809,21 @@ using LinearAlgebra ) - @timeit SimMetaData.HourGlass "05 Update To Half TimeStep" HalfTimeStep(SimMetaData, SimConstants, SimParticles, Positionₙ⁺, Velocityₙ⁺, ρₙ⁺, dρdtI, dt₂) - - - @timeit SimMetaData.HourGlass "06 Half LimitDensityAtBoundary" LimitDensityAtBoundary!(ρₙ⁺, SimConstants.ρ₀, MotionLimiter) - @timeit SimMetaData.HourGlass "Motion" ProgressMotion(SimParticles, dt₂, MotionDefinition, SimMetaData) - @timeit SimMetaData.HourGlass "07 Pressure" Pressure!(SimParticles.Pressure, ρₙ⁺,SimConstants) + @timeit SimMetaData.HourGlass "07 Pressure" PressureHalfStep!(SimParticles.Pressure, Density, dρdtI, MotionLimiter, SimConstants.ρ₀, dt₂, SimConstants) @timeit SimMetaData.HourGlass "08 Second NeighborLoop" NeighborLoopPerParticle!( SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData, SimConstants, SimParticles, ParticleRanges, CellDict, NeighborCellLists, dρdtI, Acceleration, ∇Cᵢ, ∇◌rᵢ, AccelerationMax, - Position = Positionₙ⁺, - Density = ρₙ⁺, - Velocity = Velocityₙ⁺, + Position = HalfStepPosition(Position, Velocity, MotionLimiter, dt₂), + Density = HalfStepDensity(Density, dρdtI, MotionLimiter, SimConstants.ρ₀, dt₂), + Velocity = HalfStepVelocity(Velocity, Acceleration, MotionLimiter, dt₂), ) @timeit SimMetaData.HourGlass "09 Final LimitDensityAtBoundary" LimitDensityAtBoundary!(Density, SimConstants.ρ₀, MotionLimiter) - @timeit SimMetaData.HourGlass "10 Final Density" DensityEpsi!(Density, dρdtI, ρₙ⁺, dt) + @timeit SimMetaData.HourGlass "10 Final Density" DensityEpsi!(Density, dρdtI, MotionLimiter, SimConstants.ρ₀, dt, dt₂) @timeit SimMetaData.HourGlass "11 Update To Final TimeStep" FullTimeStep(SimMetaData, SimKernel, SimConstants, SimParticles, ∇Cᵢ, ∇◌rᵢ, dt) diff --git a/src/SimulationEquations.jl b/src/SimulationEquations.jl index fead76e3..596121d1 100644 --- a/src/SimulationEquations.jl +++ b/src/SimulationEquations.jl @@ -23,6 +23,17 @@ end end end +@inline function PressureHalfStep!(Press, Density, dρdtI, MotionLimiter, ρ₀, Δt₂, SimulationConstants) + @unpack c₀ = SimulationConstants + @inbounds for i ∈ eachindex(Press, Density, dρdtI, MotionLimiter) + ρₙ⁺ = Density[i] + dρdtI[i] * Δt₂ + if (ρₙ⁺ < ρ₀) * !Bool(MotionLimiter[i]) + ρₙ⁺ = ρ₀ + end + Press[i] = EquationOfStateGamma7(ρₙ⁺, c₀, ρ₀) + end +end + # This is to handle the special factor multiplied on density in the time stepping procedure, when # using symplectic time stepping @inline function DensityEpsi!(Density, dρdtIₙ⁺,ρₙ⁺,Δt) @@ -32,6 +43,17 @@ end end end +@inline function DensityEpsi!(Density, dρdtI, MotionLimiter, ρ₀, Δt, Δt₂) + @inbounds for i in eachindex(Density, dρdtI, MotionLimiter) + ρₙ⁺ = Density[i] + dρdtI[i] * Δt₂ + if (ρₙ⁺ < ρ₀) * !Bool(MotionLimiter[i]) + ρₙ⁺ = ρ₀ + end + epsi = - (dρdtI[i] / ρₙ⁺) * Δt + Density[i] *= (2 - epsi) / (2 + epsi) + end +end + # This version of the function using !Bool(MotionLimiter) instead of BoundaryBool @inline function LimitDensityAtBoundary!(Density,ρ₀, MotionLimiter) @inbounds for i in eachindex(Density) From f372950d94445f7a335e6b8f53147cbf7f55d23c Mon Sep 17 00:00:00 2001 From: AhmedSalih3d <36305327+AhmedSalih3d@users.noreply.github.com> Date: Sat, 17 Jan 2026 23:07:43 +0100 Subject: [PATCH 6/8] =?UTF-8?q?Allow=20Update=CE=94x!=20with=20lazy=20vect?= =?UTF-8?q?ors?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- src/SPHNeighborList.jl | 15 +++++++++++++++ 1 file changed, 15 insertions(+) diff --git a/src/SPHNeighborList.jl b/src/SPHNeighborList.jl index 1e84c7fa..b57b7c65 100644 --- a/src/SPHNeighborList.jl +++ b/src/SPHNeighborList.jl @@ -153,4 +153,19 @@ Returns the new Δx. return Δx + 4 * maxd end +@inline function UpdateΔx!(Δx::T, + posₙ⁺, + pos::AbstractVector{SVector{D, T}}) where {D, T<:Real} + maxd = zero(T) + @inbounds for i in eachindex(pos) + diff = posₙ⁺[i] - pos[i] + sumsq = dot(diff, diff) + nrm = sqrt(sumsq) + if nrm > maxd + maxd = nrm + end + end + return Δx + 4 * maxd +end + end From 1115f4bca6cc6da762b7f3f45a7dd6ad5af1c12c Mon Sep 17 00:00:00 2001 From: AhmedSalih3d <36305327+AhmedSalih3d@users.noreply.github.com> Date: Sat, 17 Jan 2026 23:11:18 +0100 Subject: [PATCH 7/8] Import dot for neighbor deltas --- src/SPHNeighborList.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/src/SPHNeighborList.jl b/src/SPHNeighborList.jl index b57b7c65..de882717 100644 --- a/src/SPHNeighborList.jl +++ b/src/SPHNeighborList.jl @@ -3,6 +3,7 @@ module SPHNeighborList export ConstructStencil, ExtractCells!, UpdateNeighbors!, BuildNeighborCellLists!, ComputeCellParticleCounts, ComputeCellNeighborCounts, UpdateΔx! using StaticArrays +using LinearAlgebra: dot function ConstructStencil(V::Val{d}) where d return CartesianIndices(ntuple(_ -> -1:1, V)) From a5ab84d4e33dc9b0c2791cd92d509bbd5c31bde3 Mon Sep 17 00:00:00 2001 From: AhmedSalih3d <36305327+AhmedSalih3d@users.noreply.github.com> Date: Sat, 17 Jan 2026 23:20:04 +0100 Subject: [PATCH 8/8] Export PressureHalfStep! --- src/SimulationEquations.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/SimulationEquations.jl b/src/SimulationEquations.jl index 596121d1..49bc02f9 100644 --- a/src/SimulationEquations.jl +++ b/src/SimulationEquations.jl @@ -1,6 +1,6 @@ module SimulationEquations -export EquationOfState, EquationOfStateGamma7, Pressure!, DensityEpsi!, LimitDensityAtBoundary!, ConstructGravitySVector, InverseHydrostaticEquationOfState, Estimate7thRoot +export EquationOfState, EquationOfStateGamma7, Pressure!, PressureHalfStep!, DensityEpsi!, LimitDensityAtBoundary!, ConstructGravitySVector, InverseHydrostaticEquationOfState, Estimate7thRoot using StaticArrays using Parameters