From b7e730b07e772d62e074cc3ed562a341c7e259cd Mon Sep 17 00:00:00 2001 From: AhmedSalih3d <36305327+AhmedSalih3d@users.noreply.github.com> Date: Thu, 22 Jan 2026 02:17:58 +0100 Subject: [PATCH 01/15] Add AdvancedMDBC option --- README.md | 4 + src/PreProcess.jl | 12 ++ src/SPHCellList.jl | 193 +++++++++++++++++++++++++ src/SPHExample.jl | 2 +- src/SimulationMetaDataConfiguration.jl | 3 +- 5 files changed, 212 insertions(+), 2 deletions(-) diff --git a/README.md b/README.md index 4d4dc776..19862028 100644 --- a/README.md +++ b/README.md @@ -88,6 +88,10 @@ and adjust the simulation parameters or the `ComputerInteractions!` function. Run the script to start the simulation. Results are written in `hdfvtk` format which can be loaded with ParaView 5.12 or newer. Output is written asynchronously, so files finish flushing when the simulation completes. +Boundary conditions can be selected with the `MDBCMode` type parameter in +`SimulationMetaData`. `SimpleMDBC` remains available alongside the newer +`AdvancedMDBC` option for the updated numerical checks, pressure cloning, and +no-slip velocity treatment. To color exported cell grids by particle counts, set `ExportGridCellParticleCounts=true` in `SimulationMetaData`. This also adds a `ParticleNeighborsPerCell` array that includes each cell's particle count minus 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..be6dd5f4 100644 --- a/src/SPHCellList.jl +++ b/src/SPHCellList.jl @@ -617,6 +617,51 @@ using LinearAlgebra return bΔ, AΔ end + Base.@propagate_inbounds function ComputeInteractionsMDBCAdvanced!(SimKernel, SimConstants, + Position, Density, Velocity, + ParticleType, GhostPoints, i, j) + @unpack m₀ = SimConstants + @unpack h⁻¹, H² = SimKernel + + Dimensions = length(first(Position)) + DimensionsPlus = Dimensions + 1 + bΔ = zero(SVector{DimensionsPlus, eltype(Density)}) + AΔ = zero(SMatrix{DimensionsPlus, DimensionsPlus, eltype(Density)}) + kernel_sum = zero(eltype(Density)) + velocity_sum = zero(eltype(Velocity)) + div_pos = zero(eltype(Density)) + + if ParticleType[j] == Fluid + xᵢⱼ = GhostPoints[i] - Position[j] + xᵢⱼ² = dot(xᵢⱼ, xᵢⱼ) + if xᵢⱼ² <= H² + dᵢⱼ = sqrt(abs(xᵢⱼ²)) + q = clamp(dᵢⱼ * h⁻¹, 0.0, 2.0) + + ρⱼ = Density[j] + Wᵢⱼ = @fastpow SPHKernels.Wᵢⱼ(SimKernel, q) + ∇ᵢWᵢⱼ = @fastpow ∇Wᵢⱼ(SimKernel, q, xᵢⱼ) + + Vⱼ = m₀ / ρⱼ + VⱼWᵢⱼ = Vⱼ * Wᵢⱼ + + bΔ = SVector{DimensionsPlus, eltype(Density)}(m₀ * Wᵢⱼ, (m₀ * ∇ᵢWᵢⱼ)...) + xⱼᵢ = -xᵢⱼ + first_column = [VⱼWᵢⱼ; Vⱼ * ∇ᵢWᵢⱼ] + AΔ = SMatrix{DimensionsPlus, DimensionsPlus, eltype(Density), DimensionsPlus*DimensionsPlus}( + first_column..., + ((xⱼᵢ * first_column')')... + ) + + kernel_sum = VⱼWᵢⱼ + velocity_sum = Velocity[j] * VⱼWᵢⱼ + div_pos = Vⱼ * dot(Position[j] - GhostPoints[i], ∇ᵢWᵢⱼ) + end + end + + return bΔ, AΔ, kernel_sum, velocity_sum, div_pos + end + function ApplyMDBCBeforeHalf!(::SimulationMetaData{D,T,S,K,NoMDBC,L}, _args...) where {D,T,S<:ShiftingMode, K<:KernelOutputMode, L<:LogMode} return nothing end @@ -637,6 +682,30 @@ using LinearAlgebra return nothing 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 + ) 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)) + KernelSums = @alloc(T, length(SimParticles.Position)) + VelocitySums = @alloc(eltype(SimParticles.Velocity), length(SimParticles.Position)) + DivPos = @alloc(T, length(SimParticles.Position)) + UniqueCellsView = view(UniqueCells, 1:SimMetaData.IndexCounter) + NeighborLoopMDBCAdvanced!(SimKernel, SimConstants, ParticleRanges, UniqueCellsView, + SimParticles, bᵧ, Aᵧ, KernelSums, VelocitySums, DivPos; + Position = Position, Density = Density, Velocity = Velocity) + ApplyMDBCCorrectionAdvanced(SimConstants, SimParticles, bᵧ, Aᵧ, KernelSums, VelocitySums, DivPos) + end + + return nothing + end + function ApplyMDBCCorrection(SimConstants, SimParticles, bᵧ, Aᵧ) Position = SimParticles.Position @@ -662,6 +731,122 @@ using LinearAlgebra end end end + + function NeighborLoopMDBCAdvanced!(SimKernel, SimConstants, ParticleRanges, UniqueCellsView, + SimParticles, bᵧ, Aᵧ, KernelSums, VelocitySums, DivPos; + Position = SimParticles.Position, + Density = SimParticles.Density, + Velocity = SimParticles.Velocity) + @unpack GhostPoints = SimParticles + ParticleType = SimParticles.Type + FullStencil = ConstructStencil(Val(length(first(Position)))) + + @inbounds @threads for iter in eachindex(GhostPoints) + GhostPoint = GhostPoints[iter] + if !iszero(GhostPoint) + b_acc = zero(bᵧ[iter]) + A_acc = zero(Aᵧ[iter]) + kernel_sum_acc = zero(KernelSums[iter]) + velocity_sum_acc = zero(VelocitySums[iter]) + div_pos_acc = zero(DivPos[iter]) + + GhostCellIndex = f(SimKernel, GhostPoint) + @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Δ, kernel_sum, velocity_sum, div_pos = ComputeInteractionsMDBCAdvanced!( + SimKernel, SimConstants, Position, Density, Velocity, ParticleType, GhostPoints, iter, j + ) + b_acc += bΔ + A_acc += AΔ + kernel_sum_acc += kernel_sum + velocity_sum_acc += velocity_sum + div_pos_acc += div_pos + end + end + + bᵧ[iter] = b_acc + Aᵧ[iter] = A_acc + KernelSums[iter] = kernel_sum_acc + VelocitySums[iter] = velocity_sum_acc + DivPos[iter] = div_pos_acc + end + end + + return nothing + end + + function ApplyMDBCCorrectionAdvanced(SimConstants, SimParticles, bᵧ, Aᵧ, KernelSums, VelocitySums, DivPos) + @unpack Position, Density, GhostPoints, GhostNormals, Velocity, Acceleration, Pressure = SimParticles + @unpack ρ₀, c₀, Cb⁻¹, g = SimConstants + + kernel_sum_threshold = eltype(KernelSums)(0.1) + det_threshold = eltype(KernelSums)(1e-3) + cond_threshold = eltype(KernelSums)(50.0) + + @inbounds @simd ivdep for i in eachindex(Position) + if iszero(GhostPoints[i]) + continue + end + + if DivPos[i] <= zero(DivPos[i]) + continue + end + + A = Aᵧ[i] + kernel_sum = KernelSums[i] + ghost_density = ρ₀ + grad_density = zero(SVector{length(Position[i]), eltype(ρ₀)}) + + use_matrix = false + if kernel_sum >= kernel_sum_threshold && abs(det(A)) >= det_threshold + use_matrix = cond(A) < cond_threshold + end + + if use_matrix + ghost_state = A \ bᵧ[i] + ghost_density = first(ghost_state) + grad_density = SVector{length(ghost_state) - 1, eltype(ghost_state)}(ghost_state[2:end]) + elseif first(A) > zero(eltype(A)) + ghost_density = first(bᵧ[i]) / first(A) + end + + diff = Position[i] - GhostPoints[i] + boundary_density = ghost_density + dot(diff, grad_density) + boundary_density = isnan(boundary_density) ? ρ₀ : boundary_density + + normal = GhostNormals[i] + if !iszero(normal) + normal_distance = dot(GhostPoints[i] - Position[i], normal) + gravity_vec = ConstructGravitySVector(normal, g) + gravity_normal = dot(gravity_vec, normal) + acceleration_normal = dot(Acceleration[i], normal) + ghost_pressure = EquationOfStateGamma7(ghost_density, c₀, ρ₀) + boundary_pressure = ghost_pressure + ρ₀ * (gravity_normal - acceleration_normal) * normal_distance + density_argument = one(eltype(boundary_pressure)) + boundary_pressure * Cb⁻¹ + if density_argument > zero(density_argument) + boundary_density = ρ₀ + InverseHydrostaticEquationOfState(ρ₀, boundary_pressure, Cb⁻¹) + Pressure[i] = boundary_pressure + else + Pressure[i] = EquationOfStateGamma7(boundary_density, c₀, ρ₀) + end + else + Pressure[i] = EquationOfStateGamma7(boundary_density, c₀, ρ₀) + end + + Density[i] = isnan(boundary_density) ? ρ₀ : boundary_density + + if kernel_sum >= kernel_sum_threshold + ghost_velocity = VelocitySums[i] / kernel_sum + prescribed_velocity = Velocity[i] + Velocity[i] = (prescribed_velocity * 2) - ghost_velocity + end + end + end function GenerateMotionDetails(SimParticles, SimGeometry, Dimensions, FloatType) # Assuming group markers are sequential @@ -774,6 +959,14 @@ using LinearAlgebra @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 Apply MDBC before Full TimeStep" ApplyMDBCBeforeHalf!( + SimMetaData, SimKernel, SimConstants, SimParticles, ParticleRanges, UniqueCells; + Position = Positionₙ⁺, + Density = ρₙ⁺, + Velocity = Velocityₙ⁺, + ) + end @timeit SimMetaData.HourGlass "08 Second NeighborLoop" NeighborLoopPerParticle!( SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData, SimConstants, SimParticles, ParticleRanges, CellListIndices, 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/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 From 252e69ddf7775ce9d8fc01f0cd1a008a04b5a8a1 Mon Sep 17 00:00:00 2001 From: AhmedSalih3d <36305327+AhmedSalih3d@users.noreply.github.com> Date: Thu, 22 Jan 2026 02:21:54 +0100 Subject: [PATCH 02/15] Fix AdvancedMDBC SIMD loop --- src/SPHCellList.jl | 90 ++++++++++++++++++++++------------------------ 1 file changed, 42 insertions(+), 48 deletions(-) diff --git a/src/SPHCellList.jl b/src/SPHCellList.jl index be6dd5f4..386f54a6 100644 --- a/src/SPHCellList.jl +++ b/src/SPHCellList.jl @@ -789,61 +789,55 @@ using LinearAlgebra cond_threshold = eltype(KernelSums)(50.0) @inbounds @simd ivdep for i in eachindex(Position) - if iszero(GhostPoints[i]) - continue - end - - if DivPos[i] <= zero(DivPos[i]) - continue - end - - A = Aᵧ[i] - kernel_sum = KernelSums[i] - ghost_density = ρ₀ - grad_density = zero(SVector{length(Position[i]), eltype(ρ₀)}) - - use_matrix = false - if kernel_sum >= kernel_sum_threshold && abs(det(A)) >= det_threshold - use_matrix = cond(A) < cond_threshold - end + if !(iszero(GhostPoints[i]) || DivPos[i] <= zero(DivPos[i])) + A = Aᵧ[i] + kernel_sum = KernelSums[i] + ghost_density = ρ₀ + grad_density = zero(SVector{length(Position[i]), eltype(ρ₀)}) + + use_matrix = false + if kernel_sum >= kernel_sum_threshold && abs(det(A)) >= det_threshold + use_matrix = cond(A) < cond_threshold + end - if use_matrix - ghost_state = A \ bᵧ[i] - ghost_density = first(ghost_state) - grad_density = SVector{length(ghost_state) - 1, eltype(ghost_state)}(ghost_state[2:end]) - elseif first(A) > zero(eltype(A)) - ghost_density = first(bᵧ[i]) / first(A) - end + if use_matrix + ghost_state = A \ bᵧ[i] + ghost_density = first(ghost_state) + grad_density = SVector{length(ghost_state) - 1, eltype(ghost_state)}(ghost_state[2:end]) + elseif first(A) > zero(eltype(A)) + ghost_density = first(bᵧ[i]) / first(A) + end - diff = Position[i] - GhostPoints[i] - boundary_density = ghost_density + dot(diff, grad_density) - boundary_density = isnan(boundary_density) ? ρ₀ : boundary_density - - normal = GhostNormals[i] - if !iszero(normal) - normal_distance = dot(GhostPoints[i] - Position[i], normal) - gravity_vec = ConstructGravitySVector(normal, g) - gravity_normal = dot(gravity_vec, normal) - acceleration_normal = dot(Acceleration[i], normal) - ghost_pressure = EquationOfStateGamma7(ghost_density, c₀, ρ₀) - boundary_pressure = ghost_pressure + ρ₀ * (gravity_normal - acceleration_normal) * normal_distance - density_argument = one(eltype(boundary_pressure)) + boundary_pressure * Cb⁻¹ - if density_argument > zero(density_argument) - boundary_density = ρ₀ + InverseHydrostaticEquationOfState(ρ₀, boundary_pressure, Cb⁻¹) - Pressure[i] = boundary_pressure + diff = Position[i] - GhostPoints[i] + boundary_density = ghost_density + dot(diff, grad_density) + boundary_density = isnan(boundary_density) ? ρ₀ : boundary_density + + normal = GhostNormals[i] + if !iszero(normal) + normal_distance = dot(GhostPoints[i] - Position[i], normal) + gravity_vec = ConstructGravitySVector(normal, g) + gravity_normal = dot(gravity_vec, normal) + acceleration_normal = dot(Acceleration[i], normal) + ghost_pressure = EquationOfStateGamma7(ghost_density, c₀, ρ₀) + boundary_pressure = ghost_pressure + ρ₀ * (gravity_normal - acceleration_normal) * normal_distance + density_argument = one(eltype(boundary_pressure)) + boundary_pressure * Cb⁻¹ + if density_argument > zero(density_argument) + boundary_density = ρ₀ + InverseHydrostaticEquationOfState(ρ₀, boundary_pressure, Cb⁻¹) + Pressure[i] = boundary_pressure + else + Pressure[i] = EquationOfStateGamma7(boundary_density, c₀, ρ₀) + end else Pressure[i] = EquationOfStateGamma7(boundary_density, c₀, ρ₀) end - else - Pressure[i] = EquationOfStateGamma7(boundary_density, c₀, ρ₀) - end - Density[i] = isnan(boundary_density) ? ρ₀ : boundary_density + Density[i] = isnan(boundary_density) ? ρ₀ : boundary_density - if kernel_sum >= kernel_sum_threshold - ghost_velocity = VelocitySums[i] / kernel_sum - prescribed_velocity = Velocity[i] - Velocity[i] = (prescribed_velocity * 2) - ghost_velocity + if kernel_sum >= kernel_sum_threshold + ghost_velocity = VelocitySums[i] / kernel_sum + prescribed_velocity = Velocity[i] + Velocity[i] = (prescribed_velocity * 2) - ghost_velocity + end end end end From 74fa28397dda9ae9e04e26c236bb6421485c0384 Mon Sep 17 00:00:00 2001 From: AhmedSalih3d <36305327+AhmedSalih3d@users.noreply.github.com> Date: Thu, 22 Jan 2026 02:26:37 +0100 Subject: [PATCH 03/15] Guard AdvancedMDBC against NaNs --- src/SPHCellList.jl | 39 ++++++++++++++++++++++++--------------- 1 file changed, 24 insertions(+), 15 deletions(-) diff --git a/src/SPHCellList.jl b/src/SPHCellList.jl index 386f54a6..9d1f6c0c 100644 --- a/src/SPHCellList.jl +++ b/src/SPHCellList.jl @@ -796,34 +796,41 @@ using LinearAlgebra grad_density = zero(SVector{length(Position[i]), eltype(ρ₀)}) use_matrix = false - if kernel_sum >= kernel_sum_threshold && abs(det(A)) >= det_threshold - use_matrix = cond(A) < cond_threshold + if isfinite(kernel_sum) && kernel_sum >= kernel_sum_threshold && all(isfinite, A) && abs(det(A)) >= det_threshold + condition_number = cond(A) + use_matrix = isfinite(condition_number) && condition_number < cond_threshold end if use_matrix ghost_state = A \ bᵧ[i] - ghost_density = first(ghost_state) - grad_density = SVector{length(ghost_state) - 1, eltype(ghost_state)}(ghost_state[2:end]) - elseif first(A) > zero(eltype(A)) + if all(isfinite, ghost_state) + ghost_density = first(ghost_state) + grad_density = SVector{length(ghost_state) - 1, eltype(ghost_state)}(ghost_state[2:end]) + end + elseif isfinite(first(A)) && first(A) > zero(eltype(A)) && isfinite(first(bᵧ[i])) ghost_density = first(bᵧ[i]) / first(A) end diff = Position[i] - GhostPoints[i] boundary_density = ghost_density + dot(diff, grad_density) - boundary_density = isnan(boundary_density) ? ρ₀ : boundary_density + boundary_density = isfinite(boundary_density) ? boundary_density : ρ₀ normal = GhostNormals[i] - if !iszero(normal) + if !iszero(normal) && all(isfinite, normal) normal_distance = dot(GhostPoints[i] - Position[i], normal) gravity_vec = ConstructGravitySVector(normal, g) gravity_normal = dot(gravity_vec, normal) acceleration_normal = dot(Acceleration[i], normal) ghost_pressure = EquationOfStateGamma7(ghost_density, c₀, ρ₀) boundary_pressure = ghost_pressure + ρ₀ * (gravity_normal - acceleration_normal) * normal_distance - density_argument = one(eltype(boundary_pressure)) + boundary_pressure * Cb⁻¹ - if density_argument > zero(density_argument) - boundary_density = ρ₀ + InverseHydrostaticEquationOfState(ρ₀, boundary_pressure, Cb⁻¹) - Pressure[i] = boundary_pressure + if isfinite(boundary_pressure) + density_argument = one(eltype(boundary_pressure)) + boundary_pressure * Cb⁻¹ + if isfinite(density_argument) && density_argument > zero(density_argument) + boundary_density = ρ₀ + InverseHydrostaticEquationOfState(ρ₀, boundary_pressure, Cb⁻¹) + Pressure[i] = boundary_pressure + else + Pressure[i] = EquationOfStateGamma7(boundary_density, c₀, ρ₀) + end else Pressure[i] = EquationOfStateGamma7(boundary_density, c₀, ρ₀) end @@ -831,12 +838,14 @@ using LinearAlgebra Pressure[i] = EquationOfStateGamma7(boundary_density, c₀, ρ₀) end - Density[i] = isnan(boundary_density) ? ρ₀ : boundary_density + Density[i] = isfinite(boundary_density) ? boundary_density : ρ₀ - if kernel_sum >= kernel_sum_threshold + if isfinite(kernel_sum) && kernel_sum >= kernel_sum_threshold ghost_velocity = VelocitySums[i] / kernel_sum - prescribed_velocity = Velocity[i] - Velocity[i] = (prescribed_velocity * 2) - ghost_velocity + if all(isfinite, ghost_velocity) + prescribed_velocity = Velocity[i] + Velocity[i] = (prescribed_velocity * 2) - ghost_velocity + end end end end From aee2a381f4d588f25090dd38acb936fff63285aa Mon Sep 17 00:00:00 2001 From: AhmedSalih3d <36305327+AhmedSalih3d@users.noreply.github.com> Date: Thu, 22 Jan 2026 02:30:56 +0100 Subject: [PATCH 04/15] Stabilize AdvancedMDBC outputs --- src/SPHCellList.jl | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/SPHCellList.jl b/src/SPHCellList.jl index 9d1f6c0c..cf1b381a 100644 --- a/src/SPHCellList.jl +++ b/src/SPHCellList.jl @@ -810,6 +810,7 @@ using LinearAlgebra elseif isfinite(first(A)) && first(A) > zero(eltype(A)) && isfinite(first(bᵧ[i])) ghost_density = first(bᵧ[i]) / first(A) end + ghost_density = isfinite(ghost_density) ? ghost_density : ρ₀ diff = Position[i] - GhostPoints[i] boundary_density = ghost_density + dot(diff, grad_density) @@ -839,6 +840,9 @@ using LinearAlgebra end Density[i] = isfinite(boundary_density) ? boundary_density : ρ₀ + if !isfinite(Pressure[i]) + Pressure[i] = EquationOfStateGamma7(Density[i], c₀, ρ₀) + end if isfinite(kernel_sum) && kernel_sum >= kernel_sum_threshold ghost_velocity = VelocitySums[i] / kernel_sum From 9c6cac1c54667b36df753c14bdc3bec7117e21e4 Mon Sep 17 00:00:00 2001 From: AhmedSalih3d <36305327+AhmedSalih3d@users.noreply.github.com> Date: Thu, 22 Jan 2026 02:39:35 +0100 Subject: [PATCH 05/15] Align AdvancedMDBC checks --- src/SPHCellList.jl | 112 +++++++++++++++++++++++++-------------------- 1 file changed, 63 insertions(+), 49 deletions(-) diff --git a/src/SPHCellList.jl b/src/SPHCellList.jl index cf1b381a..d7f63198 100644 --- a/src/SPHCellList.jl +++ b/src/SPHCellList.jl @@ -788,68 +788,82 @@ using LinearAlgebra det_threshold = eltype(KernelSums)(1e-3) cond_threshold = eltype(KernelSums)(50.0) - @inbounds @simd ivdep for i in eachindex(Position) - if !(iszero(GhostPoints[i]) || DivPos[i] <= zero(DivPos[i])) - A = Aᵧ[i] - kernel_sum = KernelSums[i] - ghost_density = ρ₀ - grad_density = zero(SVector{length(Position[i]), eltype(ρ₀)}) - - use_matrix = false - if isfinite(kernel_sum) && kernel_sum >= kernel_sum_threshold && all(isfinite, A) && abs(det(A)) >= det_threshold + @inbounds for i in eachindex(Position) + if iszero(GhostPoints[i]) + continue + end + + if !(isfinite(DivPos[i]) && DivPos[i] > zero(DivPos[i])) + Density[i] = ρ₀ + Pressure[i] = EquationOfStateGamma7(ρ₀, c₀, ρ₀) + continue + end + + A = Aᵧ[i] + kernel_sum = KernelSums[i] + ghost_density = ρ₀ + grad_density = zero(SVector{length(Position[i]), eltype(ρ₀)}) + + use_shepherd = true + use_matrix = false + + if isfinite(kernel_sum) && kernel_sum >= kernel_sum_threshold && all(isfinite, A) + det_value = det(A) + if isfinite(det_value) && abs(det_value) >= det_threshold condition_number = cond(A) use_matrix = isfinite(condition_number) && condition_number < cond_threshold + use_shepherd = !use_matrix end + end - if use_matrix - ghost_state = A \ bᵧ[i] - if all(isfinite, ghost_state) - ghost_density = first(ghost_state) - grad_density = SVector{length(ghost_state) - 1, eltype(ghost_state)}(ghost_state[2:end]) - end - elseif isfinite(first(A)) && first(A) > zero(eltype(A)) && isfinite(first(bᵧ[i])) - ghost_density = first(bᵧ[i]) / first(A) + if use_matrix + ghost_state = A \ bᵧ[i] + if all(isfinite, ghost_state) + ghost_density = first(ghost_state) + grad_density = SVector{length(ghost_state) - 1, eltype(ghost_state)}(ghost_state[2:end]) end - ghost_density = isfinite(ghost_density) ? ghost_density : ρ₀ - - diff = Position[i] - GhostPoints[i] - boundary_density = ghost_density + dot(diff, grad_density) - boundary_density = isfinite(boundary_density) ? boundary_density : ρ₀ - - normal = GhostNormals[i] - if !iszero(normal) && all(isfinite, normal) - normal_distance = dot(GhostPoints[i] - Position[i], normal) - gravity_vec = ConstructGravitySVector(normal, g) - gravity_normal = dot(gravity_vec, normal) - acceleration_normal = dot(Acceleration[i], normal) - ghost_pressure = EquationOfStateGamma7(ghost_density, c₀, ρ₀) - boundary_pressure = ghost_pressure + ρ₀ * (gravity_normal - acceleration_normal) * normal_distance - if isfinite(boundary_pressure) - density_argument = one(eltype(boundary_pressure)) + boundary_pressure * Cb⁻¹ - if isfinite(density_argument) && density_argument > zero(density_argument) - boundary_density = ρ₀ + InverseHydrostaticEquationOfState(ρ₀, boundary_pressure, Cb⁻¹) - Pressure[i] = boundary_pressure - else - Pressure[i] = EquationOfStateGamma7(boundary_density, c₀, ρ₀) - end + elseif use_shepherd && isfinite(first(A)) && first(A) > zero(eltype(A)) && isfinite(first(bᵧ[i])) + ghost_density = first(bᵧ[i]) / first(A) + end + + ghost_density = isfinite(ghost_density) ? ghost_density : ρ₀ + diff = Position[i] - GhostPoints[i] + boundary_density = ghost_density + dot(diff, grad_density) + boundary_density = isfinite(boundary_density) ? boundary_density : ρ₀ + + normal = GhostNormals[i] + if !iszero(normal) && all(isfinite, normal) + normal_distance = dot(GhostPoints[i] - Position[i], normal) + gravity_vec = ConstructGravitySVector(normal, g) + gravity_normal = dot(gravity_vec, normal) + acceleration_normal = dot(Acceleration[i], normal) + ghost_pressure = EquationOfStateGamma7(ghost_density, c₀, ρ₀) + boundary_pressure = ghost_pressure + ρ₀ * (gravity_normal - acceleration_normal) * normal_distance + if isfinite(boundary_pressure) + density_argument = one(eltype(boundary_pressure)) + boundary_pressure * Cb⁻¹ + if isfinite(density_argument) && density_argument > zero(density_argument) + boundary_density = ρ₀ + InverseHydrostaticEquationOfState(ρ₀, boundary_pressure, Cb⁻¹) + Pressure[i] = boundary_pressure else Pressure[i] = EquationOfStateGamma7(boundary_density, c₀, ρ₀) end else Pressure[i] = EquationOfStateGamma7(boundary_density, c₀, ρ₀) end + else + Pressure[i] = EquationOfStateGamma7(boundary_density, c₀, ρ₀) + end - Density[i] = isfinite(boundary_density) ? boundary_density : ρ₀ - if !isfinite(Pressure[i]) - Pressure[i] = EquationOfStateGamma7(Density[i], c₀, ρ₀) - end + Density[i] = isfinite(boundary_density) ? boundary_density : ρ₀ + if !isfinite(Pressure[i]) + Pressure[i] = EquationOfStateGamma7(Density[i], c₀, ρ₀) + end - if isfinite(kernel_sum) && kernel_sum >= kernel_sum_threshold - ghost_velocity = VelocitySums[i] / kernel_sum - if all(isfinite, ghost_velocity) - prescribed_velocity = Velocity[i] - Velocity[i] = (prescribed_velocity * 2) - ghost_velocity - end + if isfinite(kernel_sum) && kernel_sum >= kernel_sum_threshold + ghost_velocity = VelocitySums[i] / kernel_sum + if all(isfinite, ghost_velocity) + prescribed_velocity = Velocity[i] + Velocity[i] = (prescribed_velocity * 2) - ghost_velocity end end end From a240629f3c6e76994ee724d8865b91b58ae129bf Mon Sep 17 00:00:00 2001 From: AhmedSalih3d <36305327+AhmedSalih3d@users.noreply.github.com> Date: Thu, 22 Jan 2026 02:45:07 +0100 Subject: [PATCH 06/15] Guard MapFloor for nonfinite --- src/SPHNeighborList.jl | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/SPHNeighborList.jl b/src/SPHNeighborList.jl index c4965635..be9cf53d 100644 --- a/src/SPHNeighborList.jl +++ b/src/SPHNeighborList.jl @@ -64,6 +64,9 @@ Extracts the cells for each particle based on their positions and the inverse cu # This is different than just doing muladd(x,InverseCutOff,0.5) because it rounds towards zero. # Consider -1.7 + 0.5, this would give -1.2 and then truncated 1, but we want -2, therefore absolute addition beforehand # We add 0.5 instead of 1, to ensure proper rounding behavior when restoring the sign for negative numbers. + if !isfinite(X) + return zero(Int) + end Int(sign(X)) * unsafe_trunc(Int, muladd(abs(X), InverseCutOff, 0.5)) end From 152b9d4d4e3cb1aef3b604c170af1b2d092bd9c6 Mon Sep 17 00:00:00 2001 From: AhmedSalih3d <36305327+AhmedSalih3d@users.noreply.github.com> Date: Thu, 22 Jan 2026 02:53:09 +0100 Subject: [PATCH 07/15] Add BoundOnOff scaling for MDBC --- src/PreProcess.jl | 2 ++ src/SPHCellList.jl | 34 ++++++++++++++++++++++------------ src/SPHNeighborList.jl | 3 --- 3 files changed, 24 insertions(+), 15 deletions(-) diff --git a/src/PreProcess.jl b/src/PreProcess.jl index 8bbfb05c..488f95ec 100644 --- a/src/PreProcess.jl +++ b/src/PreProcess.jl @@ -128,6 +128,8 @@ function AllocateDataStructures(SimGeometry::Vector{<:Geometry{Dimensions, Float GhostNormals = zeros(PositionType, NumberOfPoints) ParticleFields = merge(ParticleFields, (; GhostNormals = GhostNormals)) end + BoundOnOff = ones(PositionUnderlyingType, NumberOfPoints) + ParticleFields = merge(ParticleFields, (; BoundOnOff = BoundOnOff)) ParticleFields = merge(ParticleFields, (; ID = Idp)) SimParticles = StructArray(ParticleFields) diff --git a/src/SPHCellList.jl b/src/SPHCellList.jl index d7f63198..3f9e0bbc 100644 --- a/src/SPHCellList.jl +++ b/src/SPHCellList.jl @@ -379,6 +379,7 @@ using LinearAlgebra SV<:SPHViscosity} @unpack m₀, dx = SimConstants @unpack h⁻¹, H², h = SimKernel + @unpack BoundOnOff = SimParticles xᵢⱼ = Position[i] - Position[j] xᵢⱼ² = dot(xᵢⱼ, xᵢⱼ) @@ -395,7 +396,8 @@ using LinearAlgebra vⱼ = Velocity[j] vᵢⱼ = vᵢ - vⱼ density_symmetric_term = dot(-vᵢⱼ, ∇ᵢWᵢⱼ) - dρdt⁺ = -ρᵢ * (m₀ / ρⱼ) * density_symmetric_term + mⱼ = m₀ * BoundOnOff[j] + dρdt⁺ = -ρᵢ * (mⱼ / ρⱼ) * density_symmetric_term Dᵢ, _ = compute_density_diffusion(SimDensityDiffusion, SimKernel, SimConstants, SimParticles, xᵢⱼ, ∇ᵢWᵢⱼ, dᵢⱼ², i, j, ParticleType) @@ -405,7 +407,7 @@ using LinearAlgebra Pⱼ = Pressure[j] Pfac = (Pᵢ + Pⱼ) / (ρᵢ * ρⱼ) f_ab = tensile_correction(SimKernel, Pᵢ, ρᵢ, Pⱼ, ρⱼ, q, dx) - dvdt⁺ = -m₀ * (Pfac + f_ab) * ∇ᵢWᵢⱼ + dvdt⁺ = -mⱼ * (Pfac + f_ab) * ∇ᵢWᵢⱼ visc_term, _ = compute_viscosity(SimViscosity, SimKernel, SimConstants, SimParticles, xᵢⱼ, vᵢⱼ, ∇ᵢWᵢⱼ, dᵢⱼ², i, j) @@ -426,6 +428,7 @@ using LinearAlgebra SV<:SPHViscosity} @unpack m₀, dx = SimConstants @unpack h⁻¹, H², h = SimKernel + @unpack BoundOnOff = SimParticles xᵢⱼ = Position[i] - Position[j] xᵢⱼ² = dot(xᵢⱼ, xᵢⱼ) @@ -442,7 +445,8 @@ using LinearAlgebra vⱼ = Velocity[j] vᵢⱼ = vᵢ - vⱼ density_symmetric_term = dot(-vᵢⱼ, ∇ᵢWᵢⱼ) - dρdt⁺ = -ρᵢ * (m₀ / ρⱼ) * density_symmetric_term + mⱼ = m₀ * BoundOnOff[j] + dρdt⁺ = -ρᵢ * (mⱼ / ρⱼ) * density_symmetric_term Dᵢ, _ = compute_density_diffusion(SimDensityDiffusion, SimKernel, SimConstants, SimParticles, xᵢⱼ, ∇ᵢWᵢⱼ, dᵢⱼ², i, j, ParticleType) @@ -452,7 +456,7 @@ using LinearAlgebra Pⱼ = Pressure[j] Pfac = (Pᵢ + Pⱼ) / (ρᵢ * ρⱼ) f_ab = tensile_correction(SimKernel, Pᵢ, ρᵢ, Pⱼ, ρⱼ, q, dx) - dvdt⁺ = -m₀ * (Pfac + f_ab) * ∇ᵢWᵢⱼ + dvdt⁺ = -mⱼ * (Pfac + f_ab) * ∇ᵢWᵢⱼ visc_term, _ = compute_viscosity(SimViscosity, SimKernel, SimConstants, SimParticles, xᵢⱼ, vᵢⱼ, ∇ᵢWᵢⱼ, dᵢⱼ², i, j) @@ -475,6 +479,7 @@ using LinearAlgebra SV<:SPHViscosity} @unpack m₀, dx = SimConstants @unpack h⁻¹, H², h = SimKernel + @unpack BoundOnOff = SimParticles xᵢⱼ = Position[i] - Position[j] xᵢⱼ² = dot(xᵢⱼ, xᵢⱼ) @@ -491,7 +496,8 @@ using LinearAlgebra vⱼ = Velocity[j] vᵢⱼ = vᵢ - vⱼ density_symmetric_term = dot(-vᵢⱼ, ∇ᵢWᵢⱼ) - dρdt⁺ = -ρᵢ * (m₀ / ρⱼ) * density_symmetric_term + mⱼ = m₀ * BoundOnOff[j] + dρdt⁺ = -ρᵢ * (mⱼ / ρⱼ) * density_symmetric_term Dᵢ, _ = compute_density_diffusion(SimDensityDiffusion, SimKernel, SimConstants, SimParticles, xᵢⱼ, ∇ᵢWᵢⱼ, dᵢⱼ², i, j, ParticleType) @@ -501,7 +507,7 @@ using LinearAlgebra Pⱼ = Pressure[j] Pfac = (Pᵢ + Pⱼ) / (ρᵢ * ρⱼ) f_ab = tensile_correction(SimKernel, Pᵢ, ρᵢ, Pⱼ, ρⱼ, q, dx) - dvdt⁺ = -m₀ * (Pfac + f_ab) * ∇ᵢWᵢⱼ + dvdt⁺ = -mⱼ * (Pfac + f_ab) * ∇ᵢWᵢⱼ visc_term, _ = compute_viscosity(SimViscosity, SimKernel, SimConstants, SimParticles, xᵢⱼ, vᵢⱼ, ∇ᵢWᵢⱼ, dᵢⱼ², i, j) @@ -511,7 +517,7 @@ using LinearAlgebra MotionLimiterCondition = MotionLimiterValue(eltype(ρᵢ), ParticleType[i]) * MotionLimiterValue(eltype(ρᵢ), ParticleType[j]) shift_c_acc += (m₀ / ρᵢ) * ∇ᵢWᵢⱼ - shift_r_acc += (m₀ / ρⱼ) * dot(-xᵢⱼ, ∇ᵢWᵢⱼ) * MotionLimiterCondition + shift_r_acc += (mⱼ / ρⱼ) * dot(-xᵢⱼ, ∇ᵢWᵢⱼ) * MotionLimiterCondition end return dρdt_acc, acc_acc, kernel_acc, kernel_grad_acc, shift_c_acc, shift_r_acc @@ -529,6 +535,7 @@ using LinearAlgebra SV<:SPHViscosity} @unpack m₀, dx = SimConstants @unpack h⁻¹, H², h = SimKernel + @unpack BoundOnOff = SimParticles xᵢⱼ = Position[i] - Position[j] xᵢⱼ² = dot(xᵢⱼ, xᵢⱼ) @@ -545,7 +552,8 @@ using LinearAlgebra vⱼ = Velocity[j] vᵢⱼ = vᵢ - vⱼ density_symmetric_term = dot(-vᵢⱼ, ∇ᵢWᵢⱼ) - dρdt⁺ = -ρᵢ * (m₀ / ρⱼ) * density_symmetric_term + mⱼ = m₀ * BoundOnOff[j] + dρdt⁺ = -ρᵢ * (mⱼ / ρⱼ) * density_symmetric_term Dᵢ, _ = compute_density_diffusion(SimDensityDiffusion, SimKernel, SimConstants, SimParticles, xᵢⱼ, ∇ᵢWᵢⱼ, dᵢⱼ², i, j, ParticleType) @@ -555,7 +563,7 @@ using LinearAlgebra Pⱼ = Pressure[j] Pfac = (Pᵢ + Pⱼ) / (ρᵢ * ρⱼ) f_ab = tensile_correction(SimKernel, Pᵢ, ρᵢ, Pⱼ, ρⱼ, q, dx) - dvdt⁺ = -m₀ * (Pfac + f_ab) * ∇ᵢWᵢⱼ + dvdt⁺ = -mⱼ * (Pfac + f_ab) * ∇ᵢWᵢⱼ visc_term, _ = compute_viscosity(SimViscosity, SimKernel, SimConstants, SimParticles, xᵢⱼ, vᵢⱼ, ∇ᵢWᵢⱼ, dᵢⱼ², i, j) @@ -563,7 +571,7 @@ using LinearAlgebra MotionLimiterCondition = MotionLimiterValue(eltype(ρᵢ), ParticleType[i]) * MotionLimiterValue(eltype(ρᵢ), ParticleType[j]) shift_c_acc += (m₀ / ρᵢ) * ∇ᵢWᵢⱼ - shift_r_acc += (m₀ / ρⱼ) * dot(-xᵢⱼ, ∇ᵢWᵢⱼ) * MotionLimiterCondition + shift_r_acc += (mⱼ / ρⱼ) * dot(-xᵢⱼ, ∇ᵢWᵢⱼ) * MotionLimiterCondition end return dρdt_acc, acc_acc, shift_c_acc, shift_r_acc @@ -781,7 +789,7 @@ using LinearAlgebra end function ApplyMDBCCorrectionAdvanced(SimConstants, SimParticles, bᵧ, Aᵧ, KernelSums, VelocitySums, DivPos) - @unpack Position, Density, GhostPoints, GhostNormals, Velocity, Acceleration, Pressure = SimParticles + @unpack Position, Density, GhostPoints, GhostNormals, Velocity, Acceleration, Pressure, BoundOnOff = SimParticles @unpack ρ₀, c₀, Cb⁻¹, g = SimConstants kernel_sum_threshold = eltype(KernelSums)(0.1) @@ -794,10 +802,12 @@ using LinearAlgebra end if !(isfinite(DivPos[i]) && DivPos[i] > zero(DivPos[i])) + BoundOnOff[i] = zero(eltype(BoundOnOff)) Density[i] = ρ₀ - Pressure[i] = EquationOfStateGamma7(ρ₀, c₀, ρ₀) + Pressure[i] = zero(eltype(Pressure)) continue end + BoundOnOff[i] = one(eltype(BoundOnOff)) A = Aᵧ[i] kernel_sum = KernelSums[i] diff --git a/src/SPHNeighborList.jl b/src/SPHNeighborList.jl index be9cf53d..c4965635 100644 --- a/src/SPHNeighborList.jl +++ b/src/SPHNeighborList.jl @@ -64,9 +64,6 @@ Extracts the cells for each particle based on their positions and the inverse cu # This is different than just doing muladd(x,InverseCutOff,0.5) because it rounds towards zero. # Consider -1.7 + 0.5, this would give -1.2 and then truncated 1, but we want -2, therefore absolute addition beforehand # We add 0.5 instead of 1, to ensure proper rounding behavior when restoring the sign for negative numbers. - if !isfinite(X) - return zero(Int) - end Int(sign(X)) * unsafe_trunc(Int, muladd(abs(X), InverseCutOff, 0.5)) end From 59a4785ba97f9cb032e5e782218d20745eda0533 Mon Sep 17 00:00:00 2001 From: AhmedSalih3d <36305327+AhmedSalih3d@users.noreply.github.com> Date: Thu, 22 Jan 2026 02:58:08 +0100 Subject: [PATCH 08/15] Scale diffusion and viscosity by BoundOnOff --- src/SPHDensityDiffusionModels.jl | 12 +++++++++--- src/SPHViscosityModels.jl | 16 +++++++++++----- 2 files changed, 20 insertions(+), 8 deletions(-) diff --git a/src/SPHDensityDiffusionModels.jl b/src/SPHDensityDiffusionModels.jl index 71651654..63a42d91 100644 --- a/src/SPHDensityDiffusionModels.jl +++ b/src/SPHDensityDiffusionModels.jl @@ -69,9 +69,11 @@ struct ZeroGravityLinearDensityDiffusion <: SPHDensityDiffusion end @unpack ρ₀, m₀, c₀, δᵩ, Cb, Cb⁻¹, γ = SimConstants @unpack h, η² = SimKernel + @unpack BoundOnOff = SimParticles ρᵢ = SimParticles.Density[i] ρⱼ = SimParticles.Density[j] + mⱼ = m₀ * BoundOnOff[j] # g == 0 => skip any hydrostatic parts @@ -80,7 +82,7 @@ struct ZeroGravityLinearDensityDiffusion <: SPHDensityDiffusion end ρⱼᵢ = ρⱼ - ρᵢ ψᵢⱼ = 2 * ρⱼᵢ * (-xᵢⱼ) * invdᵢⱼ²η² - Dᵢ = δᵩ * h * c₀ * (m₀/ρⱼ) * dot(ψᵢⱼ, ∇ᵢWᵢⱼ) + Dᵢ = δᵩ * h * c₀ * (mⱼ/ρⱼ) * dot(ψᵢⱼ, ∇ᵢWᵢⱼ) Dⱼ = -Dᵢ @@ -113,11 +115,13 @@ struct LinearDensityDiffusion <: SPHDensityDiffusion end @unpack ρ₀, m₀, c₀, δᵩ, Cb, Cb⁻¹, γ, g = SimConstants @unpack h, η² = SimKernel + @unpack BoundOnOff = SimParticles Linear_ρ_factor = (1/(Cb*γ))*ρ₀ ρᵢ = SimParticles.Density[i] ρⱼ = SimParticles.Density[j] + mⱼ = m₀ * BoundOnOff[j] Pᵢⱼᴴ = ρ₀ * (-g) * -xᵢⱼ[end] ρᵢⱼᴴ = Pᵢⱼᴴ * Linear_ρ_factor @@ -130,7 +134,7 @@ struct LinearDensityDiffusion <: SPHDensityDiffusion end MotionLimiterCondition = MotionLimiterValue(eltype(ρᵢ), ParticleType[i]) * MotionLimiterValue(eltype(ρᵢ), ParticleType[j]) - Dᵢ = δᵩ * h * c₀ * (m₀/ρⱼ) * dot(ψᵢⱼ, ∇ᵢWᵢⱼ) * MotionLimiterCondition + Dᵢ = δᵩ * h * c₀ * (mⱼ/ρⱼ) * dot(ψᵢⱼ, ∇ᵢWᵢⱼ) * MotionLimiterCondition Dⱼ = -Dᵢ return Dᵢ, Dⱼ @@ -163,9 +167,11 @@ struct ComplexDensityDiffusion <: SPHDensityDiffusion end @unpack ρ₀, m₀, c₀, δᵩ, Cb, Cb⁻¹, γ, g = SimConstants @unpack h, η² = SimKernel + @unpack BoundOnOff = SimParticles ρᵢ = SimParticles.Density[i] ρⱼ = SimParticles.Density[j] + mⱼ = m₀ * BoundOnOff[j] # In theory these two equations are not completely symmetric. # In practice it is 'good' enough and saves a lot of time to @@ -182,7 +188,7 @@ struct ComplexDensityDiffusion <: SPHDensityDiffusion end MotionLimiterCondition = MotionLimiterValue(eltype(ρᵢ), ParticleType[i]) * MotionLimiterValue(eltype(ρᵢ), ParticleType[j]) - Dᵢ = δᵩ * h * c₀ * (m₀/ρⱼ) * dot(ψᵢⱼ, ∇ᵢWᵢⱼ) * MotionLimiterCondition + Dᵢ = δᵩ * h * c₀ * (mⱼ/ρⱼ) * dot(ψᵢⱼ, ∇ᵢWᵢⱼ) * MotionLimiterCondition Dⱼ = -Dᵢ return Dᵢ, Dⱼ diff --git a/src/SPHViscosityModels.jl b/src/SPHViscosityModels.jl index 318b0419..46d7bcc3 100644 --- a/src/SPHViscosityModels.jl +++ b/src/SPHViscosityModels.jl @@ -57,16 +57,18 @@ end xᵢⱼ, vᵢⱼ, ∇ᵢWᵢⱼ, d², i, j) @unpack m₀, α, c₀ = SimConstants @unpack h, η² = SimKernel + @unpack BoundOnOff = SimParticles ρᵢ = SimParticles.Density[i] ρⱼ = SimParticles.Density[j] + mⱼ = m₀ * BoundOnOff[j] v_dot_x = dot(vᵢⱼ, xᵢⱼ) if v_dot_x < 0 ρ̄ = 0.5 * (ρᵢ + ρⱼ) μᵢⱼ = h * v_dot_x / (d² + η²) - Π = -m₀ * (-α * c₀ * μᵢⱼ) / ρ̄ * ∇ᵢWᵢⱼ + Π = -mⱼ * (-α * c₀ * μᵢⱼ) / ρ̄ * ∇ᵢWᵢⱼ return Π, -Π end @@ -77,24 +79,28 @@ end @inline function compute_viscosity(::Laminar, SimKernel, SimConstants, SimParticles, xᵢⱼ, vᵢⱼ, ∇ᵢWᵢⱼ, d², i, j) @unpack m₀, ν₀ = SimConstants @unpack η² = SimKernel + @unpack BoundOnOff = SimParticles dᵢⱼ = sqrt(abs(d²)) ρᵢ = SimParticles.Density[i] ρⱼ = SimParticles.Density[j] + mⱼ = m₀ * BoundOnOff[j] - term = (4 * m₀ * ν₀ * dot(xᵢⱼ, ∇ᵢWᵢⱼ)) / ((ρᵢ + ρⱼ) + (d² + η²)) + term = (4 * mⱼ * ν₀ * dot(xᵢⱼ, ∇ᵢWᵢⱼ)) / ((ρᵢ + ρⱼ) + (d² + η²)) return term * vᵢⱼ, -term * vᵢⱼ end # LaminarSPS: with sub-grid scale stresses. @inline function compute_viscosity(::LaminarSPS, SimKernel, SimConstants, SimParticles, xᵢⱼ, vᵢⱼ, ∇ᵢWᵢⱼ, d², i, j) @unpack m₀, dx, SmagorinskyConstant, BlinConstant = SimConstants + @unpack BoundOnOff = SimParticles t1,t2 = compute_viscosity(Laminar(), SimKernel, SimConstants, SimParticles, xᵢⱼ, vᵢⱼ, ∇ᵢWᵢⱼ, d², i, j) ρᵢ = SimParticles.Density[i] ρⱼ = SimParticles.Density[j] + mⱼ = m₀ * BoundOnOff[j] vᵢ = SimParticles.Velocity[i] vⱼ = SimParticles.Velocity[j] @@ -107,19 +113,19 @@ end # 0.0 0.0 0.0 # 0.0 0.0 0.0 # Strain *rate* tensor is the gradient of velocity - Sᵢ = ∇vᵢ = (m₀/ρⱼ) * (vⱼ - vᵢ) * ∇ᵢWᵢⱼ' + Sᵢ = ∇vᵢ = (mⱼ/ρⱼ) * (vⱼ - vᵢ) * ∇ᵢWᵢⱼ' norm_Sᵢ = sqrt(2 * sum(Sᵢ .^ 2)) νtᵢ = (SmagorinskyConstant * dx)^2 * norm_Sᵢ trace_Sᵢ = sum(diag(Sᵢ)) τᶿᵢ = 2*νtᵢ*ρᵢ * (Sᵢ - (1/3) * trace_Sᵢ * Iᴹ) - (2/3) * ρᵢ * BlinConstant * dx^2 * norm_Sᵢ^2 * Iᴹ - Sⱼ = ∇vⱼ = (m₀/ρᵢ) * (vᵢ - vⱼ) * -∇ᵢWᵢⱼ' + Sⱼ = ∇vⱼ = (mⱼ/ρᵢ) * (vᵢ - vⱼ) * -∇ᵢWᵢⱼ' norm_Sⱼ = sqrt(2 * sum(Sⱼ .^ 2)) νtⱼ = (SmagorinskyConstant * dx)^2 * norm_Sⱼ trace_Sⱼ = sum(diag(Sⱼ)) τᶿⱼ = 2*νtⱼ*ρⱼ * (Sⱼ - (1/3) * trace_Sⱼ * Iᴹ) - (2/3) * ρⱼ * BlinConstant * dx^2 * norm_Sⱼ^2 * Iᴹ # MATHEMATICALLY THIS IS DOT PRODUCT TO GO FROM TENSOR TO VECTOR, BUT USE * IN JULIA TO REPRESENT IT - dτdtᵢ = (m₀/(ρⱼ * ρᵢ)) * (τᶿᵢ + τᶿⱼ) * ∇ᵢWᵢⱼ + dτdtᵢ = (mⱼ/(ρⱼ * ρᵢ)) * (τᶿᵢ + τᶿⱼ) * ∇ᵢWᵢⱼ dτdtⱼ = -dτdtᵢ return t1 + dτdtᵢ, t2 + dτdtⱼ From 7eff8f552bf151e4fa9bb5761aeba1f23b842fa7 Mon Sep 17 00:00:00 2001 From: AhmedSalih3d <36305327+AhmedSalih3d@users.noreply.github.com> Date: Thu, 22 Jan 2026 03:05:17 +0100 Subject: [PATCH 09/15] Avoid NaN checks in MDBC correction --- src/SPHCellList.jl | 34 +++++++++++++++++++--------------- 1 file changed, 19 insertions(+), 15 deletions(-) diff --git a/src/SPHCellList.jl b/src/SPHCellList.jl index 3f9e0bbc..1129dd81 100644 --- a/src/SPHCellList.jl +++ b/src/SPHCellList.jl @@ -801,7 +801,7 @@ using LinearAlgebra continue end - if !(isfinite(DivPos[i]) && DivPos[i] > zero(DivPos[i])) + if !(DivPos[i] > zero(DivPos[i])) BoundOnOff[i] = zero(eltype(BoundOnOff)) Density[i] = ρ₀ Pressure[i] = zero(eltype(Pressure)) @@ -817,41 +817,45 @@ using LinearAlgebra use_shepherd = true use_matrix = false - if isfinite(kernel_sum) && kernel_sum >= kernel_sum_threshold && all(isfinite, A) + if kernel_sum >= kernel_sum_threshold det_value = det(A) - if isfinite(det_value) && abs(det_value) >= det_threshold + if abs(det_value) >= det_threshold condition_number = cond(A) - use_matrix = isfinite(condition_number) && condition_number < cond_threshold + use_matrix = condition_number < cond_threshold use_shepherd = !use_matrix end end if use_matrix ghost_state = A \ bᵧ[i] - if all(isfinite, ghost_state) + if !any(isnan, ghost_state) ghost_density = first(ghost_state) grad_density = SVector{length(ghost_state) - 1, eltype(ghost_state)}(ghost_state[2:end]) end - elseif use_shepherd && isfinite(first(A)) && first(A) > zero(eltype(A)) && isfinite(first(bᵧ[i])) + elseif use_shepherd && first(A) > zero(eltype(A)) && !isnan(first(bᵧ[i])) ghost_density = first(bᵧ[i]) / first(A) end - ghost_density = isfinite(ghost_density) ? ghost_density : ρ₀ + if isnan(ghost_density) + ghost_density = ρ₀ + end diff = Position[i] - GhostPoints[i] boundary_density = ghost_density + dot(diff, grad_density) - boundary_density = isfinite(boundary_density) ? boundary_density : ρ₀ + if isnan(boundary_density) + boundary_density = ρ₀ + end normal = GhostNormals[i] - if !iszero(normal) && all(isfinite, normal) + if !iszero(normal) normal_distance = dot(GhostPoints[i] - Position[i], normal) gravity_vec = ConstructGravitySVector(normal, g) gravity_normal = dot(gravity_vec, normal) acceleration_normal = dot(Acceleration[i], normal) ghost_pressure = EquationOfStateGamma7(ghost_density, c₀, ρ₀) boundary_pressure = ghost_pressure + ρ₀ * (gravity_normal - acceleration_normal) * normal_distance - if isfinite(boundary_pressure) + if !isnan(boundary_pressure) density_argument = one(eltype(boundary_pressure)) + boundary_pressure * Cb⁻¹ - if isfinite(density_argument) && density_argument > zero(density_argument) + if density_argument > zero(density_argument) boundary_density = ρ₀ + InverseHydrostaticEquationOfState(ρ₀, boundary_pressure, Cb⁻¹) Pressure[i] = boundary_pressure else @@ -864,14 +868,14 @@ using LinearAlgebra Pressure[i] = EquationOfStateGamma7(boundary_density, c₀, ρ₀) end - Density[i] = isfinite(boundary_density) ? boundary_density : ρ₀ - if !isfinite(Pressure[i]) + Density[i] = isnan(boundary_density) ? ρ₀ : boundary_density + if isnan(Pressure[i]) Pressure[i] = EquationOfStateGamma7(Density[i], c₀, ρ₀) end - if isfinite(kernel_sum) && kernel_sum >= kernel_sum_threshold + if kernel_sum >= kernel_sum_threshold ghost_velocity = VelocitySums[i] / kernel_sum - if all(isfinite, ghost_velocity) + if !any(isnan, ghost_velocity) prescribed_velocity = Velocity[i] Velocity[i] = (prescribed_velocity * 2) - ghost_velocity end From 9e056fcdeec0fc09e761bbd519834aebb854367b Mon Sep 17 00:00:00 2001 From: AhmedSalih3d <36305327+AhmedSalih3d@users.noreply.github.com> Date: Thu, 22 Jan 2026 12:55:52 +0100 Subject: [PATCH 10/15] Scope BoundOnOff to boundary particles --- src/SPHCellList.jl | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/src/SPHCellList.jl b/src/SPHCellList.jl index 1129dd81..f6ffc703 100644 --- a/src/SPHCellList.jl +++ b/src/SPHCellList.jl @@ -789,7 +789,7 @@ using LinearAlgebra end function ApplyMDBCCorrectionAdvanced(SimConstants, SimParticles, bᵧ, Aᵧ, KernelSums, VelocitySums, DivPos) - @unpack Position, Density, GhostPoints, GhostNormals, Velocity, Acceleration, Pressure, BoundOnOff = SimParticles + @unpack Position, Density, GhostPoints, GhostNormals, Velocity, Acceleration, Pressure, BoundOnOff, Type = SimParticles @unpack ρ₀, c₀, Cb⁻¹, g = SimConstants kernel_sum_threshold = eltype(KernelSums)(0.1) @@ -797,14 +797,20 @@ using LinearAlgebra cond_threshold = eltype(KernelSums)(50.0) @inbounds for i in eachindex(Position) + if Type[i] == Fluid + BoundOnOff[i] = one(eltype(BoundOnOff)) + continue + end + if iszero(GhostPoints[i]) + BoundOnOff[i] = one(eltype(BoundOnOff)) continue end if !(DivPos[i] > zero(DivPos[i])) BoundOnOff[i] = zero(eltype(BoundOnOff)) Density[i] = ρ₀ - Pressure[i] = zero(eltype(Pressure)) + Pressure[i] = EquationOfStateGamma7(ρ₀, c₀, ρ₀) continue end BoundOnOff[i] = one(eltype(BoundOnOff)) From 51aea3d70144cf5e3ae59c0d1af6fc5006799b19 Mon Sep 17 00:00:00 2001 From: AhmedSalih3d <36305327+AhmedSalih3d@users.noreply.github.com> Date: Thu, 22 Jan 2026 13:07:27 +0100 Subject: [PATCH 11/15] Align AdvancedMDBC decision flow --- src/SPHCellList.jl | 27 ++++++--------------------- 1 file changed, 6 insertions(+), 21 deletions(-) diff --git a/src/SPHCellList.jl b/src/SPHCellList.jl index f6ffc703..5b5085a9 100644 --- a/src/SPHCellList.jl +++ b/src/SPHCellList.jl @@ -789,7 +789,7 @@ using LinearAlgebra end function ApplyMDBCCorrectionAdvanced(SimConstants, SimParticles, bᵧ, Aᵧ, KernelSums, VelocitySums, DivPos) - @unpack Position, Density, GhostPoints, GhostNormals, Velocity, Acceleration, Pressure, BoundOnOff, Type = SimParticles + @unpack Position, Density, GhostPoints, GhostNormals, Velocity, Acceleration, Pressure, BoundOnOff = SimParticles @unpack ρ₀, c₀, Cb⁻¹, g = SimConstants kernel_sum_threshold = eltype(KernelSums)(0.1) @@ -797,13 +797,7 @@ using LinearAlgebra cond_threshold = eltype(KernelSums)(50.0) @inbounds for i in eachindex(Position) - if Type[i] == Fluid - BoundOnOff[i] = one(eltype(BoundOnOff)) - continue - end - if iszero(GhostPoints[i]) - BoundOnOff[i] = one(eltype(BoundOnOff)) continue end @@ -820,16 +814,15 @@ using LinearAlgebra ghost_density = ρ₀ grad_density = zero(SVector{length(Position[i]), eltype(ρ₀)}) - use_shepherd = true + use_shepherd = kernel_sum < kernel_sum_threshold use_matrix = false - - if kernel_sum >= kernel_sum_threshold + if !use_shepherd det_value = det(A) if abs(det_value) >= det_threshold condition_number = cond(A) use_matrix = condition_number < cond_threshold - use_shepherd = !use_matrix end + use_shepherd = !use_matrix end if use_matrix @@ -838,8 +831,8 @@ using LinearAlgebra ghost_density = first(ghost_state) grad_density = SVector{length(ghost_state) - 1, eltype(ghost_state)}(ghost_state[2:end]) end - elseif use_shepherd && first(A) > zero(eltype(A)) && !isnan(first(bᵧ[i])) - ghost_density = first(bᵧ[i]) / first(A) + elseif use_shepherd && kernel_sum > zero(eltype(kernel_sum)) + ghost_density = first(bᵧ[i]) / kernel_sum end if isnan(ghost_density) @@ -1000,14 +993,6 @@ using LinearAlgebra @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 Apply MDBC before Full TimeStep" ApplyMDBCBeforeHalf!( - SimMetaData, SimKernel, SimConstants, SimParticles, ParticleRanges, UniqueCells; - Position = Positionₙ⁺, - Density = ρₙ⁺, - Velocity = Velocityₙ⁺, - ) - end @timeit SimMetaData.HourGlass "08 Second NeighborLoop" NeighborLoopPerParticle!( SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData, SimConstants, SimParticles, ParticleRanges, CellListIndices, From b3b091e26c22a259a5081bc7cba514b4dbbb1ac6 Mon Sep 17 00:00:00 2001 From: AhmedSalih3d <36305327+AhmedSalih3d@users.noreply.github.com> Date: Thu, 22 Jan 2026 13:38:39 +0100 Subject: [PATCH 12/15] Tidy SimpleMDBC README note --- README.md | 3 +- src/PreProcess.jl | 12 -- src/SPHCellList.jl | 185 ++++--------------------- src/SPHExample.jl | 2 +- src/SimulationMetaDataConfiguration.jl | 3 +- 5 files changed, 31 insertions(+), 174 deletions(-) diff --git a/README.md b/README.md index 19862028..7ac7d317 100644 --- a/README.md +++ b/README.md @@ -89,8 +89,7 @@ Run the script to start the simulation. Results are written in `hdfvtk` format which can be loaded with ParaView 5.12 or newer. Output is written asynchronously, so files finish flushing when the simulation completes. Boundary conditions can be selected with the `MDBCMode` type parameter in -`SimulationMetaData`. `SimpleMDBC` remains available alongside the newer -`AdvancedMDBC` option for the updated numerical checks, pressure cloning, and +`SimulationMetaData`. `SimpleMDBC` now includes updated numerical checks, pressure cloning, and no-slip velocity treatment. To color exported cell grids by particle counts, set `ExportGridCellParticleCounts=true` in `SimulationMetaData`. This also adds a diff --git a/src/PreProcess.jl b/src/PreProcess.jl index 488f95ec..4a999ad7 100644 --- a/src/PreProcess.jl +++ b/src/PreProcess.jl @@ -220,16 +220,4 @@ 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 5b5085a9..0b8f6559 100644 --- a/src/SPHCellList.jl +++ b/src/SPHCellList.jl @@ -299,7 +299,7 @@ using LinearAlgebra function NeighborLoopMDBC!(SimKernel, SimMetaData::SimulationMetaData{Dimensions, FloatType, SMode, KMode, BMode, LMode}, SimConstants, ParticleRanges, UniqueCellsView, - SimParticles, bᵧ, Aᵧ) where {Dimensions, FloatType, SMode, KMode, BMode, LMode} + SimParticles, bᵧ, Aᵧ, KernelSums, DivPos) where {Dimensions, FloatType, SMode, KMode, BMode, LMode} @unpack Position, Density, GhostPoints, GhostNormals = SimParticles ParticleType = SimParticles.Type @@ -313,6 +313,8 @@ using LinearAlgebra # zero‐initialize per‐ghost accumulators b_acc = zero(bᵧ[iter]) # an SVector{D+1,FloatType} A_acc = zero(Aᵧ[iter]) # an SMatrix{D+1,D+1,FloatType} + kernel_sum_acc = zero(KernelSums[iter]) + div_pos_acc = zero(DivPos[iter]) # compute and accumulate into the locals GhostCellIndex = f(SimKernel, GhostPoints[iter]) @@ -329,17 +331,21 @@ using LinearAlgebra for j in StartIndex_:EndIndex_ # change ComputeInteractions to take & return contributions, e.g.: - bΔ, AΔ = ComputeInteractionsMDBC!(SimKernel, SimMetaData, SimConstants, - Position, Density, ParticleType, - GhostPoints, iter, j) + bΔ, AΔ, kernel_sum, div_pos = ComputeInteractionsMDBC!(SimKernel, SimMetaData, SimConstants, + Position, Density, ParticleType, + GhostPoints, iter, j) b_acc += bΔ A_acc += AΔ + kernel_sum_acc += kernel_sum + div_pos_acc += div_pos end end # write out once bᵧ[iter] = b_acc Aᵧ[iter] = A_acc + KernelSums[iter] = kernel_sum_acc + DivPos[iter] = div_pos_acc end end @@ -586,6 +592,8 @@ using LinearAlgebra # always zero‐initialize bΔ = zero(SVector{DimensionsPlus,FloatType}) AΔ = zero(SMatrix{DimensionsPlus, DimensionsPlus,FloatType}) + kernel_sum = zero(FloatType) + div_pos = zero(FloatType) # ᵢ is ghost node! ⱼ is fluid node @@ -618,56 +626,14 @@ using LinearAlgebra first_column..., ((xⱼᵢ * first_column')')... ) - end - end - - - return bΔ, AΔ - end - - Base.@propagate_inbounds function ComputeInteractionsMDBCAdvanced!(SimKernel, SimConstants, - Position, Density, Velocity, - ParticleType, GhostPoints, i, j) - @unpack m₀ = SimConstants - @unpack h⁻¹, H² = SimKernel - - Dimensions = length(first(Position)) - DimensionsPlus = Dimensions + 1 - bΔ = zero(SVector{DimensionsPlus, eltype(Density)}) - AΔ = zero(SMatrix{DimensionsPlus, DimensionsPlus, eltype(Density)}) - kernel_sum = zero(eltype(Density)) - velocity_sum = zero(eltype(Velocity)) - div_pos = zero(eltype(Density)) - - if ParticleType[j] == Fluid - xᵢⱼ = GhostPoints[i] - Position[j] - xᵢⱼ² = dot(xᵢⱼ, xᵢⱼ) - if xᵢⱼ² <= H² - dᵢⱼ = sqrt(abs(xᵢⱼ²)) - q = clamp(dᵢⱼ * h⁻¹, 0.0, 2.0) - - ρⱼ = Density[j] - Wᵢⱼ = @fastpow SPHKernels.Wᵢⱼ(SimKernel, q) - ∇ᵢWᵢⱼ = @fastpow ∇Wᵢⱼ(SimKernel, q, xᵢⱼ) - - Vⱼ = m₀ / ρⱼ - VⱼWᵢⱼ = Vⱼ * Wᵢⱼ - - bΔ = SVector{DimensionsPlus, eltype(Density)}(m₀ * Wᵢⱼ, (m₀ * ∇ᵢWᵢⱼ)...) - xⱼᵢ = -xᵢⱼ - first_column = [VⱼWᵢⱼ; Vⱼ * ∇ᵢWᵢⱼ] - AΔ = SMatrix{DimensionsPlus, DimensionsPlus, eltype(Density), DimensionsPlus*DimensionsPlus}( - first_column..., - ((xⱼᵢ * first_column')')... - ) kernel_sum = VⱼWᵢⱼ - velocity_sum = Velocity[j] * VⱼWᵢⱼ div_pos = Vⱼ * dot(Position[j] - GhostPoints[i], ∇ᵢWᵢⱼ) end end - - return bΔ, AΔ, kernel_sum, velocity_sum, div_pos + + + return bΔ, AΔ, kernel_sum, div_pos end function ApplyMDBCBeforeHalf!(::SimulationMetaData{D,T,S,K,NoMDBC,L}, _args...) where {D,T,S<:ShiftingMode, K<:KernelOutputMode, L<:LogMode} @@ -678,138 +644,51 @@ using LinearAlgebra SimKernel, SimConstants, SimParticles, ParticleRanges, UniqueCells ) 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)) - UniqueCellsView = view(UniqueCells, 1:SimMetaData.IndexCounter) - NeighborLoopMDBC!(SimKernel, SimMetaData, SimConstants, ParticleRanges, UniqueCellsView, SimParticles, bᵧ, Aᵧ) - ApplyMDBCCorrection(SimConstants, SimParticles, bᵧ, Aᵧ) - end - - return nothing - 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 - ) 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)) KernelSums = @alloc(T, length(SimParticles.Position)) - VelocitySums = @alloc(eltype(SimParticles.Velocity), length(SimParticles.Position)) DivPos = @alloc(T, length(SimParticles.Position)) UniqueCellsView = view(UniqueCells, 1:SimMetaData.IndexCounter) - NeighborLoopMDBCAdvanced!(SimKernel, SimConstants, ParticleRanges, UniqueCellsView, - SimParticles, bᵧ, Aᵧ, KernelSums, VelocitySums, DivPos; - Position = Position, Density = Density, Velocity = Velocity) - ApplyMDBCCorrectionAdvanced(SimConstants, SimParticles, bᵧ, Aᵧ, KernelSums, VelocitySums, DivPos) + NeighborLoopMDBC!(SimKernel, SimMetaData, SimConstants, ParticleRanges, UniqueCellsView, SimParticles, bᵧ, Aᵧ, KernelSums, DivPos) + ApplyMDBCCorrection(SimConstants, SimParticles, bᵧ, Aᵧ, KernelSums, DivPos) end return nothing end - - function ApplyMDBCCorrection(SimConstants, SimParticles, bᵧ, Aᵧ) + function ApplyMDBCCorrection(SimConstants, SimParticles, bᵧ, Aᵧ, KernelSums, DivPos) Position = SimParticles.Position Density = SimParticles.Density GhostPoints = SimParticles.GhostPoints + BoundOnOff = SimParticles.BoundOnOff + Pressure = SimParticles.Pressure + GhostNormals = SimParticles.GhostNormals + Acceleration = SimParticles.Acceleration ρ₀ = SimConstants.ρ₀ - #https://github.com/DualSPHysics/DualSPHysics/blob/f4fa76ad5083873fa1c6dd3b26cdce89c55a9aeb/src/source/JSphCpu_mdbc.cpp#L347 - @inbounds @simd ivdep for i in eachindex(Position) - A = Aᵧ[i] - - # Since Aᵧ is not reset anymore, we need to check if it is zero - if !iszero(GhostPoints[i]) - if abs(det(A)) >= 1e-3 - GhostPointDensity = A \ bᵧ[i] - diff = Position[i] - GhostPoints[i] - v1 = first(GhostPointDensity) + sum(GhostPointDensity[j+1] * diff[j] for j in eachindex(diff)) - Density[i] = isnan(v1) ? ρ₀ : v1 - elseif first(A) > 0.0 - v = first(bᵧ[i]) / first(A) - Density[i] = isnan(v) ? ρ₀ : v - end - end - end - end - - function NeighborLoopMDBCAdvanced!(SimKernel, SimConstants, ParticleRanges, UniqueCellsView, - SimParticles, bᵧ, Aᵧ, KernelSums, VelocitySums, DivPos; - Position = SimParticles.Position, - Density = SimParticles.Density, - Velocity = SimParticles.Velocity) - @unpack GhostPoints = SimParticles - ParticleType = SimParticles.Type - FullStencil = ConstructStencil(Val(length(first(Position)))) - - @inbounds @threads for iter in eachindex(GhostPoints) - GhostPoint = GhostPoints[iter] - if !iszero(GhostPoint) - b_acc = zero(bᵧ[iter]) - A_acc = zero(Aᵧ[iter]) - kernel_sum_acc = zero(KernelSums[iter]) - velocity_sum_acc = zero(VelocitySums[iter]) - div_pos_acc = zero(DivPos[iter]) - - GhostCellIndex = f(SimKernel, GhostPoint) - @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Δ, kernel_sum, velocity_sum, div_pos = ComputeInteractionsMDBCAdvanced!( - SimKernel, SimConstants, Position, Density, Velocity, ParticleType, GhostPoints, iter, j - ) - b_acc += bΔ - A_acc += AΔ - kernel_sum_acc += kernel_sum - velocity_sum_acc += velocity_sum - div_pos_acc += div_pos - end - end - - bᵧ[iter] = b_acc - Aᵧ[iter] = A_acc - KernelSums[iter] = kernel_sum_acc - VelocitySums[iter] = velocity_sum_acc - DivPos[iter] = div_pos_acc - end - end - - return nothing - end - - function ApplyMDBCCorrectionAdvanced(SimConstants, SimParticles, bᵧ, Aᵧ, KernelSums, VelocitySums, DivPos) - @unpack Position, Density, GhostPoints, GhostNormals, Velocity, Acceleration, Pressure, BoundOnOff = SimParticles - @unpack ρ₀, c₀, Cb⁻¹, g = SimConstants + c₀ = SimConstants.c₀ + Cb⁻¹ = SimConstants.Cb⁻¹ + g = SimConstants.g kernel_sum_threshold = eltype(KernelSums)(0.1) det_threshold = eltype(KernelSums)(1e-3) cond_threshold = eltype(KernelSums)(50.0) - @inbounds for i in eachindex(Position) + @inbounds @simd ivdep for i in eachindex(Position) + A = Aᵧ[i] + if iszero(GhostPoints[i]) continue end if !(DivPos[i] > zero(DivPos[i])) BoundOnOff[i] = zero(eltype(BoundOnOff)) - Density[i] = ρ₀ - Pressure[i] = EquationOfStateGamma7(ρ₀, c₀, ρ₀) continue end BoundOnOff[i] = one(eltype(BoundOnOff)) - A = Aᵧ[i] kernel_sum = KernelSums[i] ghost_density = ρ₀ grad_density = zero(SVector{length(Position[i]), eltype(ρ₀)}) @@ -871,14 +750,6 @@ using LinearAlgebra if isnan(Pressure[i]) Pressure[i] = EquationOfStateGamma7(Density[i], c₀, ρ₀) end - - if kernel_sum >= kernel_sum_threshold - ghost_velocity = VelocitySums[i] / kernel_sum - if !any(isnan, ghost_velocity) - prescribed_velocity = Velocity[i] - Velocity[i] = (prescribed_velocity * 2) - ghost_velocity - end - end end end diff --git a/src/SPHExample.jl b/src/SPHExample.jl index d687b29e..7929bdf3 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, AdvancedMDBC, + MDBCMode, NoMDBC, SimpleMDBC, LogMode, NoLog, StoreLog, TimeSteppingMode, SymplecticTimeStepping, SingleNeighborTimeStepping diff --git a/src/SimulationMetaDataConfiguration.jl b/src/SimulationMetaDataConfiguration.jl index 6b783f0e..b92381bf 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, AdvancedMDBC, + MDBCMode, NoMDBC, SimpleMDBC, LogMode, NoLog, StoreLog, TimeSteppingMode, SymplecticTimeStepping, SingleNeighborTimeStepping @@ -21,7 +21,6 @@ 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 From 8fcd9b9aaec3772f36c4e61f48befac3e7cacdc0 Mon Sep 17 00:00:00 2001 From: AhmedSalih3d <36305327+AhmedSalih3d@users.noreply.github.com> Date: Thu, 22 Jan 2026 13:42:17 +0100 Subject: [PATCH 13/15] Drop simd from MDBC correction --- src/SPHCellList.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/SPHCellList.jl b/src/SPHCellList.jl index 0b8f6559..09fc93cc 100644 --- a/src/SPHCellList.jl +++ b/src/SPHCellList.jl @@ -676,7 +676,7 @@ using LinearAlgebra det_threshold = eltype(KernelSums)(1e-3) cond_threshold = eltype(KernelSums)(50.0) - @inbounds @simd ivdep for i in eachindex(Position) + @inbounds for i in eachindex(Position) A = Aᵧ[i] if iszero(GhostPoints[i]) From 8dd927423599e42347319e00b33e3f863b19a469 Mon Sep 17 00:00:00 2001 From: AhmedSalih3d <36305327+AhmedSalih3d@users.noreply.github.com> Date: Thu, 22 Jan 2026 13:47:37 +0100 Subject: [PATCH 14/15] Reset boundary state when MDBC off --- src/SPHCellList.jl | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/SPHCellList.jl b/src/SPHCellList.jl index 09fc93cc..2e2756ae 100644 --- a/src/SPHCellList.jl +++ b/src/SPHCellList.jl @@ -685,6 +685,8 @@ using LinearAlgebra if !(DivPos[i] > zero(DivPos[i])) BoundOnOff[i] = zero(eltype(BoundOnOff)) + Density[i] = ρ₀ + Pressure[i] = EquationOfStateGamma7(ρ₀, c₀, ρ₀) continue end BoundOnOff[i] = one(eltype(BoundOnOff)) From f073c44628bbbded83bad6e33d1d384444c5a0d1 Mon Sep 17 00:00:00 2001 From: AhmedSalih3d <36305327+AhmedSalih3d@users.noreply.github.com> Date: Thu, 22 Jan 2026 13:51:00 +0100 Subject: [PATCH 15/15] Limit BoundOnOff scaling to boundaries --- src/SPHCellList.jl | 8 ++++---- src/SPHDensityDiffusionModels.jl | 6 +++--- src/SPHViscosityModels.jl | 6 +++--- 3 files changed, 10 insertions(+), 10 deletions(-) diff --git a/src/SPHCellList.jl b/src/SPHCellList.jl index 2e2756ae..ab5712db 100644 --- a/src/SPHCellList.jl +++ b/src/SPHCellList.jl @@ -402,7 +402,7 @@ using LinearAlgebra vⱼ = Velocity[j] vᵢⱼ = vᵢ - vⱼ density_symmetric_term = dot(-vᵢⱼ, ∇ᵢWᵢⱼ) - mⱼ = m₀ * BoundOnOff[j] + mⱼ = ParticleType[j] == Fluid ? m₀ : m₀ * BoundOnOff[j] dρdt⁺ = -ρᵢ * (mⱼ / ρⱼ) * density_symmetric_term Dᵢ, _ = compute_density_diffusion(SimDensityDiffusion, SimKernel, SimConstants, SimParticles, xᵢⱼ, ∇ᵢWᵢⱼ, dᵢⱼ², i, j, ParticleType) @@ -451,7 +451,7 @@ using LinearAlgebra vⱼ = Velocity[j] vᵢⱼ = vᵢ - vⱼ density_symmetric_term = dot(-vᵢⱼ, ∇ᵢWᵢⱼ) - mⱼ = m₀ * BoundOnOff[j] + mⱼ = ParticleType[j] == Fluid ? m₀ : m₀ * BoundOnOff[j] dρdt⁺ = -ρᵢ * (mⱼ / ρⱼ) * density_symmetric_term Dᵢ, _ = compute_density_diffusion(SimDensityDiffusion, SimKernel, SimConstants, SimParticles, xᵢⱼ, ∇ᵢWᵢⱼ, dᵢⱼ², i, j, ParticleType) @@ -502,7 +502,7 @@ using LinearAlgebra vⱼ = Velocity[j] vᵢⱼ = vᵢ - vⱼ density_symmetric_term = dot(-vᵢⱼ, ∇ᵢWᵢⱼ) - mⱼ = m₀ * BoundOnOff[j] + mⱼ = ParticleType[j] == Fluid ? m₀ : m₀ * BoundOnOff[j] dρdt⁺ = -ρᵢ * (mⱼ / ρⱼ) * density_symmetric_term Dᵢ, _ = compute_density_diffusion(SimDensityDiffusion, SimKernel, SimConstants, SimParticles, xᵢⱼ, ∇ᵢWᵢⱼ, dᵢⱼ², i, j, ParticleType) @@ -558,7 +558,7 @@ using LinearAlgebra vⱼ = Velocity[j] vᵢⱼ = vᵢ - vⱼ density_symmetric_term = dot(-vᵢⱼ, ∇ᵢWᵢⱼ) - mⱼ = m₀ * BoundOnOff[j] + mⱼ = ParticleType[j] == Fluid ? m₀ : m₀ * BoundOnOff[j] dρdt⁺ = -ρᵢ * (mⱼ / ρⱼ) * density_symmetric_term Dᵢ, _ = compute_density_diffusion(SimDensityDiffusion, SimKernel, SimConstants, SimParticles, xᵢⱼ, ∇ᵢWᵢⱼ, dᵢⱼ², i, j, ParticleType) diff --git a/src/SPHDensityDiffusionModels.jl b/src/SPHDensityDiffusionModels.jl index 63a42d91..d0c1ba8f 100644 --- a/src/SPHDensityDiffusionModels.jl +++ b/src/SPHDensityDiffusionModels.jl @@ -73,7 +73,7 @@ struct ZeroGravityLinearDensityDiffusion <: SPHDensityDiffusion end ρᵢ = SimParticles.Density[i] ρⱼ = SimParticles.Density[j] - mⱼ = m₀ * BoundOnOff[j] + mⱼ = ParticleType[j] == Fluid ? m₀ : m₀ * BoundOnOff[j] # g == 0 => skip any hydrostatic parts @@ -121,7 +121,7 @@ struct LinearDensityDiffusion <: SPHDensityDiffusion end ρᵢ = SimParticles.Density[i] ρⱼ = SimParticles.Density[j] - mⱼ = m₀ * BoundOnOff[j] + mⱼ = ParticleType[j] == Fluid ? m₀ : m₀ * BoundOnOff[j] Pᵢⱼᴴ = ρ₀ * (-g) * -xᵢⱼ[end] ρᵢⱼᴴ = Pᵢⱼᴴ * Linear_ρ_factor @@ -171,7 +171,7 @@ struct ComplexDensityDiffusion <: SPHDensityDiffusion end ρᵢ = SimParticles.Density[i] ρⱼ = SimParticles.Density[j] - mⱼ = m₀ * BoundOnOff[j] + mⱼ = ParticleType[j] == Fluid ? m₀ : m₀ * BoundOnOff[j] # In theory these two equations are not completely symmetric. # In practice it is 'good' enough and saves a lot of time to diff --git a/src/SPHViscosityModels.jl b/src/SPHViscosityModels.jl index 46d7bcc3..37d02097 100644 --- a/src/SPHViscosityModels.jl +++ b/src/SPHViscosityModels.jl @@ -61,7 +61,7 @@ end ρᵢ = SimParticles.Density[i] ρⱼ = SimParticles.Density[j] - mⱼ = m₀ * BoundOnOff[j] + mⱼ = SimParticles.Type[j] == Fluid ? m₀ : m₀ * BoundOnOff[j] v_dot_x = dot(vᵢⱼ, xᵢⱼ) if v_dot_x < 0 @@ -84,7 +84,7 @@ end dᵢⱼ = sqrt(abs(d²)) ρᵢ = SimParticles.Density[i] ρⱼ = SimParticles.Density[j] - mⱼ = m₀ * BoundOnOff[j] + mⱼ = SimParticles.Type[j] == Fluid ? m₀ : m₀ * BoundOnOff[j] term = (4 * mⱼ * ν₀ * dot(xᵢⱼ, ∇ᵢWᵢⱼ)) / ((ρᵢ + ρⱼ) + (d² + η²)) return term * vᵢⱼ, -term * vᵢⱼ @@ -100,7 +100,7 @@ end ρᵢ = SimParticles.Density[i] ρⱼ = SimParticles.Density[j] - mⱼ = m₀ * BoundOnOff[j] + mⱼ = SimParticles.Type[j] == Fluid ? m₀ : m₀ * BoundOnOff[j] vᵢ = SimParticles.Velocity[i] vⱼ = SimParticles.Velocity[j]