From c642ad1cbb788d83bd494313c6a6fcd83e62f6b5 Mon Sep 17 00:00:00 2001 From: AhmedSalih3d <36305327+AhmedSalih3d@users.noreply.github.com> Date: Thu, 22 Jan 2026 13:44:36 +0100 Subject: [PATCH] Add AdvancedMDBC option --- README.md | 4 +- src/PreProcess.jl | 12 ++ src/SPHCellList.jl | 169 ++++++++++++++++++++++++- src/SPHExample.jl | 2 +- src/SimulationEquations.jl | 13 ++ src/SimulationMetaDataConfiguration.jl | 3 +- 6 files changed, 197 insertions(+), 6 deletions(-) diff --git a/README.md b/README.md index 4d4dc776..af417b31 100644 --- a/README.md +++ b/README.md @@ -21,7 +21,9 @@ The project demonstrates how to assemble a small SPH solver with Julia. It focus - **Multi-threaded execution** – achieved by spawning the neighbour loop. - **Per-particle compute loops** – per-particle loops handle threading without chunk metadata. -- **Dynamic boundary condition** – inspired by DualSPHysics. +- **Dynamic boundary condition** – inspired by DualSPHysics with `SimpleMDBC()` or + `AdvancedMDBC()` options (advanced adds numerical checks, pressure cloning, + and no-slip velocity handling). - **Density diffusion** – based on Fourtakas et al. 2019 to reduce pressure noise. - **Wendland quintic kernel** – simple and stable without tensile corrections. - **Symplectic time stepping** – choose between symplectic two-loop and single-loop midpoint updates. diff --git a/src/PreProcess.jl b/src/PreProcess.jl index 1e7630ca..8bbfb05c 100644 --- a/src/PreProcess.jl +++ b/src/PreProcess.jl @@ -218,4 +218,16 @@ function LoadMDBCNormals!(::SimulationMetaData{D,T,S,K,SimpleMDBC,L}, SimParticl return nothing end +function LoadMDBCNormals!(::SimulationMetaData{D,T,S,K,AdvancedMDBC,L}, SimParticles, path) where {D,T,S<:ShiftingMode, K<:KernelOutputMode, L<:LogMode} + if isnothing(path) + return nothing + end + _, GhostPoints, GhostNormals = LoadBoundaryNormals(Val(D), T, path) + for gi ∈ eachindex(GhostPoints) + SimParticles.GhostPoints[gi] = GhostPoints[gi] + SimParticles.GhostNormals[gi] = GhostNormals[gi] + end + return nothing +end + end diff --git a/src/SPHCellList.jl b/src/SPHCellList.jl index 1cb20127..a818501a 100644 --- a/src/SPHCellList.jl +++ b/src/SPHCellList.jl @@ -346,6 +346,55 @@ using LinearAlgebra return nothing end + function NeighborLoopMDBCAdvanced!(SimKernel, + SimMetaData::SimulationMetaData{Dimensions, FloatType, SMode, KMode, AdvancedMDBC, LMode}, + SimConstants, ParticleRanges, UniqueCellsView, + SimParticles, bᵧ, Aᵧ, VelocitySum; + Position = SimParticles.Position, + Density = SimParticles.Density, + Velocity = SimParticles.Velocity) where {Dimensions, FloatType, SMode, KMode, LMode} + + GhostPoints = SimParticles.GhostPoints + ParticleType = SimParticles.Type + + FullStencil = ConstructStencil(Val(Dimensions)) + + @inbounds @threads for iter in eachindex(GhostPoints) + GhostPoint = GhostPoints[iter] + + if !iszero(GhostPoint) + b_acc = zero(bᵧ[iter]) + A_acc = zero(Aᵧ[iter]) + VelocityAcc = zero(VelocitySum[iter]) + + GhostCellIndex = f(SimKernel, GhostPoints[iter]) + @inbounds for offset ∈ FullStencil + SCellIndex = GhostCellIndex + offset + + NeighborIdx = FindCellIndex(UniqueCellsView, SCellIndex) + + StartIndex_ = ParticleRanges[NeighborIdx] + EndIndex_ = ParticleRanges[NeighborIdx + 1] - 1 + + for j in StartIndex_:EndIndex_ + bΔ, AΔ = ComputeInteractionsMDBC!(SimKernel, SimMetaData, SimConstants, + Position, Density, ParticleType, + GhostPoints, iter, j) + b_acc += bΔ + A_acc += AΔ + VelocityAcc += AΔ[1, 1] * Velocity[j] + end + end + + bᵧ[iter] = b_acc + Aᵧ[iter] = A_acc + VelocitySum[iter] = VelocityAcc + end + end + + return nothing + end + # The previous generic `ComputeInteractions!` implementation was unused # in favour of the per-particle variants (ComputeInteractionsPerParticle! etc.). # It has been removed to reduce code size and avoid dead code. @@ -662,6 +711,100 @@ using LinearAlgebra end end end + + @inline function MDBCPositionDivergence(A) + div_pos = zero(eltype(A)) + @inbounds for idx in 2:size(A, 1) + div_pos += A[idx, idx] + end + return div_pos + end + + function ApplyMDBCBeforeHalf!(SimMetaData::SimulationMetaData{D,T,S,K,AdvancedMDBC,L}, + SimKernel, SimConstants, SimParticles, + ParticleRanges, UniqueCells; + Position = SimParticles.Position, + Density = SimParticles.Density, + Velocity = SimParticles.Velocity, + Pressure = SimParticles.Pressure, + Acceleration = SimParticles.Acceleration + ) where {D,T,S<:ShiftingMode,K<:KernelOutputMode,L<:LogMode} + @no_escape begin + DimensionsPlus = D + 1 + bᵧ = @alloc(SVector{DimensionsPlus, T}, length(SimParticles.Position)) + Aᵧ = @alloc(SMatrix{DimensionsPlus, DimensionsPlus, T, DimensionsPlus*DimensionsPlus}, length(SimParticles.Position)) + VelocitySum = @alloc(SVector{D, T}, length(SimParticles.Position)) + UniqueCellsView = view(UniqueCells, 1:SimMetaData.IndexCounter) + NeighborLoopMDBCAdvanced!(SimKernel, SimMetaData, SimConstants, ParticleRanges, UniqueCellsView, SimParticles, bᵧ, Aᵧ, VelocitySum; + Position = Position, Density = Density, Velocity = Velocity) + ApplyMDBCCorrectionAdvanced(SimConstants, SimParticles, Position, Density, Velocity, Pressure, Acceleration, bᵧ, Aᵧ, VelocitySum) + end + + return nothing + end + + function ApplyMDBCCorrectionAdvanced(SimConstants, SimParticles, Position, Density, Velocity, Pressure, Acceleration, bᵧ, Aᵧ, VelocitySum) + + GhostPoints = SimParticles.GhostPoints + GhostNormals = SimParticles.GhostNormals + + ρ₀ = SimConstants.ρ₀ + c₀ = SimConstants.c₀ + GravityVector = ConstructGravitySVector(first(Position), SimConstants.g) + + @inbounds @simd ivdep for i in eachindex(Position) + A = Aᵧ[i] + + if !iszero(GhostPoints[i]) + KernelSum = A[1, 1] + if KernelSum <= 0 + continue + end + + DivPosition = MDBCPositionDivergence(A) + if DivPosition <= 0 + continue + end + + UseShepherd = KernelSum < 0.1 + if !UseShepherd + if abs(det(A)) <= 1e-3 + UseShepherd = true + elseif cond(A) >= 50 + UseShepherd = true + end + end + + if UseShepherd + GhostDensity = first(bᵧ[i]) / KernelSum + else + GhostPointDensity = A \ bᵧ[i] + diff = Position[i] - GhostPoints[i] + GhostDensity = first(GhostPointDensity) + sum(GhostPointDensity[j+1] * diff[j] for j in eachindex(diff)) + end + + if isnan(GhostDensity) + GhostDensity = ρ₀ + end + + P_g = c₀^2 * (GhostDensity - ρ₀) + normal = GhostNormals[i] + GhostOffset = GhostPoints[i] - Position[i] + BoundaryAcceleration = Acceleration[i] + PressureTerm = ρ₀ * dot(GravityVector - BoundaryAcceleration, normal) * dot(GhostOffset, normal) + BoundaryPressure = P_g + PressureTerm + BoundaryDensity = ρ₀ + BoundaryPressure / c₀^2 + + Density[i] = BoundaryDensity + Pressure[i] = BoundaryPressure + + GhostVelocity = VelocitySum[i] / KernelSum + BoundaryVelocity = 2 * Velocity[i] - GhostVelocity + BoundaryVelocity = BoundaryVelocity - dot(BoundaryVelocity, normal) * normal + Velocity[i] = BoundaryVelocity + end + end + end function GenerateMotionDetails(SimParticles, SimGeometry, Dimensions, FloatType) # Assuming group markers are sequential @@ -769,11 +912,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.ρ₀, ParticleType) + @timeit SimMetaData.HourGlass "06 Half LimitDensityAtBoundary" LimitDensityAtBoundary!(ρₙ⁺, SimConstants.ρ₀, ParticleType, SimMetaData) @timeit SimMetaData.HourGlass "Motion" ProgressMotion(SimParticles, dt₂, MotionDefinition, SimMetaData) @timeit SimMetaData.HourGlass "07 Pressure" Pressure!(SimParticles.Pressure, ρₙ⁺, SimConstants) + if SimMetaData isa SimulationMetaData{Dimensions, FloatType, SMode, KMode, AdvancedMDBC, LMode} + @timeit SimMetaData.HourGlass "07a Half Apply Advanced MDBC" ApplyMDBCBeforeHalf!( + SimMetaData, SimKernel, SimConstants, SimParticles, ParticleRanges, UniqueCells; + Position = Positionₙ⁺, + Density = ρₙ⁺, + Velocity = Velocityₙ⁺, + Pressure = SimParticles.Pressure, + Acceleration = SimParticles.Acceleration, + ) + end @timeit SimMetaData.HourGlass "08 Second NeighborLoop" NeighborLoopPerParticle!( SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData, SimConstants, SimParticles, ParticleRanges, CellListIndices, @@ -787,11 +940,21 @@ using LinearAlgebra @timeit SimMetaData.HourGlass "03 Update To Half TimeStep" HalfTimeStep(SimMetaData, SimConstants, SimParticles, Positionₙ⁺, Velocityₙ⁺, ρₙ⁺, dρdtI, dt₂) - @timeit SimMetaData.HourGlass "04 Half LimitDensityAtBoundary" LimitDensityAtBoundary!(ρₙ⁺, SimConstants.ρ₀, ParticleType) + @timeit SimMetaData.HourGlass "04 Half LimitDensityAtBoundary" LimitDensityAtBoundary!(ρₙ⁺, SimConstants.ρ₀, ParticleType, SimMetaData) @timeit SimMetaData.HourGlass "Motion" ProgressMotion(SimParticles, dt₂, MotionDefinition, SimMetaData) @timeit SimMetaData.HourGlass "05 Pressure" Pressure!(SimParticles.Pressure, ρₙ⁺, SimConstants) + if SimMetaData isa SimulationMetaData{Dimensions, FloatType, SMode, KMode, AdvancedMDBC, LMode} + @timeit SimMetaData.HourGlass "05a Half Apply Advanced MDBC" ApplyMDBCBeforeHalf!( + SimMetaData, SimKernel, SimConstants, SimParticles, ParticleRanges, UniqueCells; + Position = Positionₙ⁺, + Density = ρₙ⁺, + Velocity = Velocityₙ⁺, + Pressure = SimParticles.Pressure, + Acceleration = SimParticles.Acceleration, + ) + end @timeit SimMetaData.HourGlass "06 NeighborLoop" NeighborLoopPerParticle!( SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData, SimConstants, SimParticles, ParticleRanges, CellListIndices, @@ -804,7 +967,7 @@ using LinearAlgebra @timeit SimMetaData.HourGlass "07 Final Density" DensityEpsi!(SimParticles.Density, dρdtI, ρₙ⁺, dt) - @timeit SimMetaData.HourGlass "08 Final LimitDensityAtBoundary" LimitDensityAtBoundary!(SimParticles.Density, SimConstants.ρ₀, ParticleType) + @timeit SimMetaData.HourGlass "08 Final LimitDensityAtBoundary" LimitDensityAtBoundary!(SimParticles.Density, SimConstants.ρ₀, ParticleType, SimMetaData) @timeit SimMetaData.HourGlass "09 Update To Final TimeStep" FullTimeStep(SimMetaData, SimKernel, SimConstants, SimParticles, Velocityₙ⁺, ∇Cᵢ, ∇◌rᵢ, dt) diff --git a/src/SPHExample.jl b/src/SPHExample.jl index 7929bdf3..d687b29e 100644 --- a/src/SPHExample.jl +++ b/src/SPHExample.jl @@ -50,7 +50,7 @@ module SPHExample using .SimulationMetaDataConfiguration export SimulationMetaData, ShiftingMode, NoShifting, PlanarShifting, KernelOutputMode, NoKernelOutput, StoreKernelOutput, - MDBCMode, NoMDBC, SimpleMDBC, + MDBCMode, NoMDBC, SimpleMDBC, AdvancedMDBC, LogMode, NoLog, StoreLog, TimeSteppingMode, SymplecticTimeStepping, SingleNeighborTimeStepping diff --git a/src/SimulationEquations.jl b/src/SimulationEquations.jl index 44b026d5..3465b868 100644 --- a/src/SimulationEquations.jl +++ b/src/SimulationEquations.jl @@ -6,6 +6,7 @@ using StaticArrays using Parameters using FastPow using ..SimulationGeometry +using ..SimulationMetaDataConfiguration @inline function EquationOfStateGamma7(ρ,c₀,ρ₀) return @fastpow ((c₀^2*ρ₀)/7) * ((ρ/ρ₀)^7 - 1) @@ -42,6 +43,18 @@ end end end +@inline function LimitDensityAtBoundary!(Density, ρ₀, ParticleType, ::MDBCMode) + return LimitDensityAtBoundary!(Density, ρ₀, ParticleType) +end + +@inline function LimitDensityAtBoundary!(_Density, _ρ₀, _ParticleType, ::AdvancedMDBC) + return nothing +end + +@inline function LimitDensityAtBoundary!(Density, ρ₀, ParticleType, SimMetaData::SimulationMetaData{D,T,S,K,B,L}) where {D,T,S,K,B<:MDBCMode,L} + return LimitDensityAtBoundary!(Density, ρ₀, ParticleType, B()) +end + @inline function ConstructGravitySVector(_::SVector{N, T}, value) where {N, T} return SVector{N, T}(ntuple(i -> i == N ? value : 0, N)) end diff --git a/src/SimulationMetaDataConfiguration.jl b/src/SimulationMetaDataConfiguration.jl index b92381bf..6b783f0e 100644 --- a/src/SimulationMetaDataConfiguration.jl +++ b/src/SimulationMetaDataConfiguration.jl @@ -6,7 +6,7 @@ using TimerOutputs export SimulationMetaData, UpdateMetaData!, ShiftingMode, NoShifting, PlanarShifting, KernelOutputMode, NoKernelOutput, StoreKernelOutput, - MDBCMode, NoMDBC, SimpleMDBC, + MDBCMode, NoMDBC, SimpleMDBC, AdvancedMDBC, LogMode, NoLog, StoreLog, TimeSteppingMode, SymplecticTimeStepping, SingleNeighborTimeStepping @@ -21,6 +21,7 @@ struct StoreKernelOutput <: KernelOutputMode end abstract type MDBCMode end struct NoMDBC <: MDBCMode end struct SimpleMDBC <: MDBCMode end +struct AdvancedMDBC <: MDBCMode end abstract type LogMode end struct NoLog <: LogMode end