From 703813414cbb89d86ac6da86142774a814a35564 Mon Sep 17 00:00:00 2001 From: AhmedSalih3d <36305327+AhmedSalih3d@users.noreply.github.com> Date: Sat, 17 Jan 2026 23:24:20 +0100 Subject: [PATCH 1/2] Make StillWedge example portable --- example/StillWedgeMDBC.jl | 11 ++--- src/SPHCellList.jl | 84 +++++++++++++++++++++++++++++--------- src/SPHNeighborList.jl | 22 +++++++++- src/SimulationEquations.jl | 24 ++++++++++- 4 files changed, 112 insertions(+), 29 deletions(-) diff --git a/example/StillWedgeMDBC.jl b/example/StillWedgeMDBC.jl index 76795e98..8e6d7168 100644 --- a/example/StillWedgeMDBC.jl +++ b/example/StillWedgeMDBC.jl @@ -7,9 +7,10 @@ let SimConstantsWedge = SimulationConstants{FloatType}(dx=0.02,c₀=42.48576250492629, δᵩ = 0.1, CFL=0.5) # SimConstantsWedge = SimulationConstants{FloatType}(dx=0.01,c₀=43.4, δᵩ = 0.1, CFL=0.2) # + OutputDirectory = joinpath(pwd(), "output", "StillWedge2D_MDBC") SimMetaDataWedge = SimulationMetaData{Dimensions,FloatType,NoShifting,NoKernelOutput,SimpleMDBC,StoreLog}( SimulationName="StillWedge", - SaveLocation="W:/Simulations/StillWedge2D_MDBC", + SaveLocation=OutputDirectory, SimulationTime=4.0, OutputTimes=0.01, VisualizeInParaview=true, @@ -19,9 +20,7 @@ let ) # If save directory is not already made, make it - if !isdir(SimMetaDataWedge.SaveLocation) - mkdir(SimMetaDataWedge.SaveLocation) - end + mkpath(SimMetaDataWedge.SaveLocation) # Assuming SimConstantsWedge is defined somewhere else with the field `dx` FixedBoundary = Geometry{Dimensions, FloatType}( @@ -49,7 +48,7 @@ let SimKernel = SPHKernelInstance{Dimensions, FloatType}(WendlandC2(); dx = SimConstantsWedge.dx) - @profview RunSimulation( + RunSimulation( SimGeometry = SimulationGeometry, SimMetaData = SimMetaDataWedge, SimConstants = SimConstantsWedge, @@ -97,5 +96,3 @@ let # display(plt) end - - diff --git a/src/SPHCellList.jl b/src/SPHCellList.jl index eb55a820..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, @@ -688,9 +731,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{ @@ -717,14 +758,25 @@ using LinearAlgebra dt = SimConstants.CFL * (SimKernel.h / SimConstants.c₀) @no_escape begin - AccelerationMax = @alloc(FloatType, length(Position)) + PositionType = eltype(Position) + PositionUnderlyingType = eltype(PositionType) + AccelerationMax = similar(Density) + dρdtI = 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) + fill!(∇Cᵢ, zero(PositionType)) + fill!(∇◌rᵢ, zero(PositionUnderlyingType)) + end while SimMetaData.TotalTime <= next_output_time(SimMetaData) @timeit SimMetaData.HourGlass "01 Calculate IndexCounter" begin - SimMetaData.Δx = UpdateΔx!(SimMetaData.Δx, Positionₙ⁺, SimParticles.Position) - ShouldRebuild = SimMetaData.Δx >= SimKernel.h + 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) @@ -757,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) @@ -805,8 +852,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 +904,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) diff --git a/src/SPHNeighborList.jl b/src/SPHNeighborList.jl index 48d245ea..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)) @@ -77,6 +78,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 +102,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 @@ -149,4 +154,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 diff --git a/src/SimulationEquations.jl b/src/SimulationEquations.jl index fead76e3..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 @@ -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 0f52844ff0da04458fa84f6d6f35aee555a6b895 Mon Sep 17 00:00:00 2001 From: AhmedSalih3d <36305327+AhmedSalih3d@users.noreply.github.com> Date: Sat, 17 Jan 2026 23:26:23 +0100 Subject: [PATCH 2/2] Add keys for half-step views --- src/SPHCellList.jl | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/SPHCellList.jl b/src/SPHCellList.jl index 95d6da91..68c1f6e9 100644 --- a/src/SPHCellList.jl +++ b/src/SPHCellList.jl @@ -42,6 +42,7 @@ using LinearAlgebra Base.IndexStyle(::Type{<:HalfStepPosition}) = IndexLinear() Base.axes(half::HalfStepPosition) = axes(half.Position) Base.size(half::HalfStepPosition) = size(half.Position) + Base.keys(half::HalfStepPosition) = keys(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} @@ -54,6 +55,7 @@ using LinearAlgebra Base.IndexStyle(::Type{<:HalfStepVelocity}) = IndexLinear() Base.axes(half::HalfStepVelocity) = axes(half.Velocity) Base.size(half::HalfStepVelocity) = size(half.Velocity) + Base.keys(half::HalfStepVelocity) = keys(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} @@ -67,6 +69,7 @@ using LinearAlgebra Base.IndexStyle(::Type{<:HalfStepDensity}) = IndexLinear() Base.axes(half::HalfStepDensity) = axes(half.Density) Base.size(half::HalfStepDensity) = size(half.Density) + Base.keys(half::HalfStepDensity) = keys(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])