From d48dcc3f98ca7db4a78f1869a1473df9db7d759a Mon Sep 17 00:00:00 2001 From: AhmedSalih3d <36305327+AhmedSalih3d@users.noreply.github.com> Date: Wed, 8 Jul 2026 21:45:07 +0200 Subject: [PATCH 1/2] Decouple output from solver state --- src/SPHCellList.jl | 54 ++++++++++++++++++++++++++++++--------------- src/TimeStepping.jl | 11 ++++++++- 2 files changed, 46 insertions(+), 19 deletions(-) diff --git a/src/SPHCellList.jl b/src/SPHCellList.jl index 3439367b..07908045 100644 --- a/src/SPHCellList.jl +++ b/src/SPHCellList.jl @@ -689,9 +689,17 @@ using TimerOutputs: @timeit end - # Per-particle local Δx removed: use single scalar `SimMetaData.Δx`. + # Runtime state that must persist across output intervals. Keeping these + # values outside SimulationLoop prevents output scheduling from changing + # neighbor-list rebuilds or adaptive time-step sequencing. + mutable struct SimulationRuntimeState{T} + Δx::T + dt::T + IsInitialized::Bool + end @inbounds function SimulationLoop(SimDensityDiffusion::SDD, SimViscosity::SV, SimKernel, + Runtime::SimulationRuntimeState{FloatType}, SimMetaData::SimulationMetaData{Dimensions, FloatType, SMode, KMode, BMode, LMode}, SimConstants, SimParticles, FullStencil, ParticleRanges, UniqueCells, CellListIndices, @@ -718,34 +726,36 @@ using TimerOutputs: @timeit ### UniqueCellsView = view(UniqueCells, 1:SimMetaData.IndexCounter) - # This code here is to initialize the first time step for each simulation loop - dt = SimConstants.CFL * (SimKernel.h / SimConstants.c₀) + dt = Runtime.dt TimeSteppingMode = SimMetaData.TimeSteppingMode @no_escape begin AccelerationMax = @alloc(FloatType, length(SimParticles.Position)) dt₂ = dt * 0.5 - SimMetaData.IndexCounter = UpdateNeighbors!(SimParticles, SimKernel.H⁻¹, SortingScratchSpace, ParticleRanges, UniqueCells, CellListIndices) - UniqueCellsView = view(UniqueCells, 1:SimMetaData.IndexCounter) - BuildNeighborCellLists!(NeighborCellLists, FullStencil, UniqueCellsView, ParticleRanges) + if !Runtime.IsInitialized + SimMetaData.IndexCounter = UpdateNeighbors!(SimParticles, SimKernel.H⁻¹, SortingScratchSpace, ParticleRanges, UniqueCells, CellListIndices) + UniqueCellsView = view(UniqueCells, 1:SimMetaData.IndexCounter) + BuildNeighborCellLists!(NeighborCellLists, FullStencil, UniqueCellsView, ParticleRanges) + Runtime.IsInitialized = true - 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 "00b Init NeighborLoop" NeighborLoopPerParticle!( - SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData, - SimConstants, SimParticles, ParticleRanges, CellListIndices, - NeighborCellLists, dρdtI, SimParticles.Acceleration, ∇Cᵢ, ∇◌rᵢ, AccelerationMax, - ) + 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 "00b Init NeighborLoop" NeighborLoopPerParticle!( + SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData, + SimConstants, SimParticles, ParticleRanges, CellListIndices, + NeighborCellLists, dρdtI, SimParticles.Acceleration, ∇Cᵢ, ∇◌rᵢ, AccelerationMax, + ) + end end NextOutputTime = next_output_time(SimMetaData) while SimMetaData.TotalTime <= NextOutputTime @timeit SimMetaData.HourGlass "01 Calculate IndexCounter" begin - SimMetaData.Δx = UpdateΔx!(SimMetaData.Δx, Positionₙ⁺, SimParticles.Position) - ShouldRebuild = SimMetaData.Δx >= SimKernel.h + Runtime.Δx = UpdateΔx!(Runtime.Δx, Positionₙ⁺, SimParticles.Position) + ShouldRebuild = Runtime.Δx >= SimKernel.h # println("Δx: ", Δx, "h: ", SimKernel.h," dt: ", SimMetaData.CurrentTimeStep, " Iteration: ", SimMetaData.Iteration, " TotalTime: ", SimMetaData.TotalTime, " OutputIterationCounter: ", SimMetaData.OutputIterationCounter) @@ -757,7 +767,7 @@ using TimerOutputs: @timeit # 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, CellListIndices) - SimMetaData.Δx = zero(eltype(dρdtI)) + Runtime.Δx = zero(eltype(dρdtI)) UniqueCellsView = view(UniqueCells, 1:SimMetaData.IndexCounter) BuildNeighborCellLists!(NeighborCellLists, FullStencil, UniqueCellsView, ParticleRanges) end @@ -819,7 +829,10 @@ using TimerOutputs: @timeit @timeit SimMetaData.HourGlass "10 Update MetaData" UpdateMetaData!(SimMetaData, dt) @timeit SimMetaData.HourGlass "11 Update TimeStep" dt = UpdateTimeStep(AccelerationMax, SimConstants, SimKernel) + dt₂ = dt * 0.5 end + + Runtime.dt = dt end return nothing @@ -913,6 +926,11 @@ using TimerOutputs: @timeit FullStencil = ConstructStencil(Val(Dimensions)) NeighborCellLists = [Int[] for _ in 1:length(UniqueCells)] _, SortingScratchSpace = Base.Sort.make_scratch(nothing, eltype(SimParticles), NumberOfPoints) + Runtime = SimulationRuntimeState( + one(FloatType) + SimKernel.h, + SimConstants.CFL * (SimKernel.h / SimConstants.c₀), + false, + ) output = SetupVTKOutput(SimMetaData, SimParticles, SimKernel, Dimensions) @@ -949,7 +967,7 @@ using TimerOutputs: @timeit @inbounds while true @timeit SimMetaData.HourGlass "00 SimulationLoop" SimulationLoop( - SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData, + SimDensityDiffusion, SimViscosity, SimKernel, Runtime, SimMetaData, SimConstants, SimParticles, FullStencil, ParticleRanges, UniqueCells, CellListIndices, SortingScratchSpace, NeighborCellLists, dρdtI, Velocityₙ⁺, Positionₙ⁺, ρₙ⁺, diff --git a/src/TimeStepping.jl b/src/TimeStepping.jl index 7973253b..9d4804af 100644 --- a/src/TimeStepping.jl +++ b/src/TimeStepping.jl @@ -33,6 +33,15 @@ function Δt(max_acceleration, SimulationConstants, SPHKernel) return CFL * min(dt_speed, dt_force) end +function Δt(_Position, _Velocity, Acceleration::AbstractVector, SimulationConstants, SPHKernel) + max_acceleration = zero(real(eltype(eltype(Acceleration)))) + @inbounds for a in Acceleration + acc_norm = norm(a) + max_acceleration = ifelse(acc_norm > max_acceleration, acc_norm, max_acceleration) + end + return Δt(max_acceleration, SimulationConstants, SPHKernel) +end + """ UpdateTimeStep(AccelerationMax, SimConstants, SimKernel) @@ -57,7 +66,7 @@ end @inline function next_output_time(times::AbstractVector, SimMetaData) idx = SimMetaData.OutputIterationCounter - if idx < length(times) + if idx <= length(times) return times[idx] else return SimMetaData.SimulationTime From d119d331da4bd46419e47ae54013ddd761efeb8d Mon Sep 17 00:00:00 2001 From: AhmedSalih3d <36305327+AhmedSalih3d@users.noreply.github.com> Date: Wed, 8 Jul 2026 22:29:54 +0200 Subject: [PATCH 2/2] Stabilize neighbor rebuild state --- src/SPHCellList.jl | 15 ++++++++++++++- src/TimeStepping.jl | 4 ++++ 2 files changed, 18 insertions(+), 1 deletion(-) diff --git a/src/SPHCellList.jl b/src/SPHCellList.jl index 07908045..06f523a4 100644 --- a/src/SPHCellList.jl +++ b/src/SPHCellList.jl @@ -737,6 +737,8 @@ using TimerOutputs: @timeit SimMetaData.IndexCounter = UpdateNeighbors!(SimParticles, SimKernel.H⁻¹, SortingScratchSpace, ParticleRanges, UniqueCells, CellListIndices) UniqueCellsView = view(UniqueCells, 1:SimMetaData.IndexCounter) BuildNeighborCellLists!(NeighborCellLists, FullStencil, UniqueCellsView, ParticleRanges) + copyto!(Positionₙ⁺, SimParticles.Position) + Runtime.Δx = zero(Runtime.Δx) Runtime.IsInitialized = true if TimeSteppingMode isa SingleNeighborTimeStepping @@ -770,6 +772,17 @@ using TimerOutputs: @timeit Runtime.Δx = zero(eltype(dρdtI)) UniqueCellsView = view(UniqueCells, 1:SimMetaData.IndexCounter) BuildNeighborCellLists!(NeighborCellLists, FullStencil, UniqueCellsView, ParticleRanges) + copyto!(Positionₙ⁺, SimParticles.Position) + + if TimeSteppingMode isa SingleNeighborTimeStepping + @timeit SimMetaData.HourGlass "01b Rebuild Pressure" Pressure!(SimParticles.Pressure, SimParticles.Density, SimConstants) + @timeit SimMetaData.HourGlass "01c Rebuild MDBC" ApplyMDBCBeforeHalf!(SimMetaData, SimKernel, SimConstants, SimParticles, ParticleRanges, UniqueCells) + @timeit SimMetaData.HourGlass "01d Rebuild NeighborLoop" NeighborLoopPerParticle!( + SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData, + SimConstants, SimParticles, ParticleRanges, CellListIndices, + NeighborCellLists, dρdtI, SimParticles.Acceleration, ∇Cᵢ, ∇◌rᵢ, AccelerationMax, + ) + end end end @@ -927,7 +940,7 @@ using TimerOutputs: @timeit NeighborCellLists = [Int[] for _ in 1:length(UniqueCells)] _, SortingScratchSpace = Base.Sort.make_scratch(nothing, eltype(SimParticles), NumberOfPoints) Runtime = SimulationRuntimeState( - one(FloatType) + SimKernel.h, + zero(FloatType), SimConstants.CFL * (SimKernel.h / SimConstants.c₀), false, ) diff --git a/src/TimeStepping.jl b/src/TimeStepping.jl index 9d4804af..b2e9d3c7 100644 --- a/src/TimeStepping.jl +++ b/src/TimeStepping.jl @@ -122,6 +122,10 @@ function HalfTimeStep(::SimulationMetaData{Dimensions, FloatType, SMode, KMode, return nothing end +function FullTimeStep(SimMetaData::SimulationMetaData, SimKernel, SimConstants, SimParticles, ∇Cᵢ, ∇◌rᵢ, dt) + return FullTimeStep(SimMetaData, SimKernel, SimConstants, SimParticles, SimParticles.Velocity, ∇Cᵢ, ∇◌rᵢ, dt) +end + function FullTimeStep(::SimulationMetaData{D,T,NoShifting,K,B,L}, SimKernel, SimConstants, SimParticles, Velocityₙ⁺, ∇Cᵢ, ∇◌rᵢ, dt) where {D,T, K<:KernelOutputMode,