From 548e212587bd5d43ba91288c4aec43009e045539 Mon Sep 17 00:00:00 2001 From: AhmedSalih3d <36305327+AhmedSalih3d@users.noreply.github.com> Date: Sun, 25 Jan 2026 19:38:53 +0100 Subject: [PATCH 1/4] Add WCSPH shifting corrections --- src/PreProcess.jl | 8 ++- src/SPHCellList.jl | 80 ++++++++++++++++---------- src/SPHExample.jl | 1 + src/SimulationMetaDataConfiguration.jl | 23 ++++++++ src/TimeStepping.jl | 46 ++++++++++++--- test/runtests.jl | 2 +- 6 files changed, 119 insertions(+), 41 deletions(-) diff --git a/src/PreProcess.jl b/src/PreProcess.jl index d00bc336..ae98da7d 100644 --- a/src/PreProcess.jl +++ b/src/PreProcess.jl @@ -178,8 +178,10 @@ function AllocateSupportDataStructures(::SimulationMetaData{D,T,NoShifting,K,B,L ∇Cᵢ = Vector{PositionType}(undef, 0) ∇◌rᵢ = Vector{PositionUnderlyingType}(undef, 0) + ∇ρᵢ = Vector{PositionType}(undef, 0) + ∇uᵢ = Vector{SMatrix{D,D,PositionUnderlyingType}}(undef, 0) - return dρdtI, Velocityₙ⁺, Positionₙ⁺, ρₙ⁺, ∇Cᵢ, ∇◌rᵢ + return dρdtI, Velocityₙ⁺, Positionₙ⁺, ρₙ⁺, ∇Cᵢ, ∇◌rᵢ, ∇ρᵢ, ∇uᵢ end function AllocateSupportDataStructures(::SimulationMetaData{D,T,S,K,B,L}, Position) where {D,T,S<:ShiftingMode, @@ -198,8 +200,10 @@ function AllocateSupportDataStructures(::SimulationMetaData{D,T,S,K,B,L}, Positi ∇Cᵢ = zeros(PositionType, NumberOfPoints) ∇◌rᵢ = zeros(PositionUnderlyingType, NumberOfPoints) + ∇ρᵢ = zeros(PositionType, NumberOfPoints) + ∇uᵢ = zeros(SMatrix{D,D,PositionUnderlyingType}, NumberOfPoints) - return dρdtI, Velocityₙ⁺, Positionₙ⁺, ρₙ⁺, ∇Cᵢ, ∇◌rᵢ + return dρdtI, Velocityₙ⁺, Positionₙ⁺, ρₙ⁺, ∇Cᵢ, ∇◌rᵢ, ∇ρᵢ, ∇uᵢ end function LoadBoundaryNormals(::Val{D}, ::Type{T}, path_mdbc) where {D, T} diff --git a/src/SPHCellList.jl b/src/SPHCellList.jl index 1cb20127..ff362d6b 100644 --- a/src/SPHCellList.jl +++ b/src/SPHCellList.jl @@ -37,7 +37,7 @@ using LinearAlgebra SimConstants, SimParticles, ParticleRanges, CellListIndices, NeighborCellLists, dρdtI, Acceleration, ∇Cᵢ, - ∇◌rᵢ, AccelerationMax; + ∇◌rᵢ, ∇ρᵢ, ∇uᵢ, AccelerationMax; Position = SimParticles.Position, Density = SimParticles.Density, Pressure = SimParticles.Pressure, @@ -93,7 +93,7 @@ using LinearAlgebra SimConstants, SimParticles, ParticleRanges, CellListIndices, NeighborCellLists, dρdtI, Acceleration, ∇Cᵢ, - ∇◌rᵢ, AccelerationMax; + ∇◌rᵢ, ∇ρᵢ, ∇uᵢ, AccelerationMax; Position = SimParticles.Position, Density = SimParticles.Density, Pressure = SimParticles.Pressure, @@ -161,7 +161,7 @@ using LinearAlgebra SimConstants, SimParticles, ParticleRanges, CellListIndices, NeighborCellLists, dρdtI, Acceleration, ∇Cᵢ, - ∇◌rᵢ, AccelerationMax; + ∇◌rᵢ, ∇ρᵢ, ∇uᵢ, AccelerationMax; Position = SimParticles.Position, Density = SimParticles.Density, Pressure = SimParticles.Pressure, @@ -175,39 +175,41 @@ using LinearAlgebra acc_acc = zero(Acceleration[i]) shift_c_acc = zero(∇Cᵢ[i]) shift_r_acc = zero(∇◌rᵢ[i]) + density_grad_acc = zero(∇ρᵢ[i]) + velocity_grad_acc = zero(∇uᵢ[i]) CellListIndex = CellListIndices[i] SameCellStart = ParticleRanges[CellListIndex] SameCellEnd = ParticleRanges[CellListIndex + 1] - 1 NeighborCellIndices = NeighborCellLists[CellListIndex] @inbounds for j in SameCellStart:(i - 1) - dρdt_acc, acc_acc, shift_c_acc, shift_r_acc = + dρdt_acc, acc_acc, shift_c_acc, shift_r_acc, density_grad_acc, velocity_grad_acc = ComputeInteractionsPerParticleNoKernel!( SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData, SimConstants, SimParticles, Position, Density, Pressure, Velocity, ParticleType, dρdt_acc, acc_acc, shift_c_acc, - shift_r_acc, i, j, + shift_r_acc, density_grad_acc, velocity_grad_acc, i, j, ) end @inbounds for j in (i + 1):SameCellEnd - dρdt_acc, acc_acc, shift_c_acc, shift_r_acc = + dρdt_acc, acc_acc, shift_c_acc, shift_r_acc, density_grad_acc, velocity_grad_acc = ComputeInteractionsPerParticleNoKernel!( SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData, SimConstants, SimParticles, Position, Density, Pressure, Velocity, ParticleType, dρdt_acc, acc_acc, shift_c_acc, - shift_r_acc, i, j, + shift_r_acc, density_grad_acc, velocity_grad_acc, i, j, ) end for NeighborIdx in NeighborCellIndices StartIndex_ = ParticleRanges[NeighborIdx] EndIndex_ = ParticleRanges[NeighborIdx + 1] - 1 @inbounds for j in StartIndex_:EndIndex_ - dρdt_acc, acc_acc, shift_c_acc, shift_r_acc = + dρdt_acc, acc_acc, shift_c_acc, shift_r_acc, density_grad_acc, velocity_grad_acc = ComputeInteractionsPerParticleNoKernel!( SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData, SimConstants, SimParticles, Position, Density, Pressure, Velocity, ParticleType, dρdt_acc, acc_acc, shift_c_acc, - shift_r_acc, i, j, + shift_r_acc, density_grad_acc, velocity_grad_acc, i, j, ) end end @@ -216,6 +218,8 @@ using LinearAlgebra Acceleration[i] = acc_acc ∇Cᵢ[i] = shift_c_acc ∇◌rᵢ[i] = shift_r_acc + ∇ρᵢ[i] = density_grad_acc + ∇uᵢ[i] = velocity_grad_acc AccelerationMax[i] = norm(acc_acc) end @@ -227,7 +231,7 @@ using LinearAlgebra SimConstants, SimParticles, ParticleRanges, CellListIndices, NeighborCellLists, dρdtI, Acceleration, ∇Cᵢ, - ∇◌rᵢ, AccelerationMax; + ∇◌rᵢ, ∇ρᵢ, ∇uᵢ, AccelerationMax; Position = SimParticles.Position, Density = SimParticles.Density, Pressure = SimParticles.Pressure, @@ -246,6 +250,8 @@ using LinearAlgebra kernel_grad_acc = zero(KernelGradient[i]) shift_c_acc = zero(∇Cᵢ[i]) shift_r_acc = zero(∇◌rᵢ[i]) + density_grad_acc = zero(∇ρᵢ[i]) + velocity_grad_acc = zero(∇uᵢ[i]) CellListIndex = CellListIndices[i] SameCellStart = ParticleRanges[CellListIndex] SameCellEnd = ParticleRanges[CellListIndex + 1] - 1 @@ -253,20 +259,22 @@ using LinearAlgebra @inbounds for j in SameCellStart:(i - 1) dρdt_acc, acc_acc, kernel_acc, kernel_grad_acc, shift_c_acc, - shift_r_acc = ComputeInteractionsPerParticle!( + shift_r_acc, density_grad_acc, velocity_grad_acc = ComputeInteractionsPerParticle!( SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData, SimConstants, SimParticles, Position, Density, Pressure, Velocity, ParticleType, dρdt_acc, acc_acc, kernel_acc, - kernel_grad_acc, shift_c_acc, shift_r_acc, i, j, + kernel_grad_acc, shift_c_acc, shift_r_acc, density_grad_acc, + velocity_grad_acc, i, j, ) end @inbounds for j in (i + 1):SameCellEnd dρdt_acc, acc_acc, kernel_acc, kernel_grad_acc, shift_c_acc, - shift_r_acc = ComputeInteractionsPerParticle!( + shift_r_acc, density_grad_acc, velocity_grad_acc = ComputeInteractionsPerParticle!( SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData, SimConstants, SimParticles, Position, Density, Pressure, Velocity, ParticleType, dρdt_acc, acc_acc, kernel_acc, - kernel_grad_acc, shift_c_acc, shift_r_acc, i, j, + kernel_grad_acc, shift_c_acc, shift_r_acc, density_grad_acc, + velocity_grad_acc, i, j, ) end for NeighborIdx in NeighborCellIndices @@ -274,11 +282,12 @@ using LinearAlgebra EndIndex_ = ParticleRanges[NeighborIdx + 1] - 1 @inbounds for j in StartIndex_:EndIndex_ dρdt_acc, acc_acc, kernel_acc, kernel_grad_acc, shift_c_acc, - shift_r_acc = ComputeInteractionsPerParticle!( + shift_r_acc, density_grad_acc, velocity_grad_acc = ComputeInteractionsPerParticle!( SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData, SimConstants, SimParticles, Position, Density, Pressure, Velocity, ParticleType, dρdt_acc, acc_acc, kernel_acc, - kernel_grad_acc, shift_c_acc, shift_r_acc, i, j, + kernel_grad_acc, shift_c_acc, shift_r_acc, density_grad_acc, + velocity_grad_acc, i, j, ) end end @@ -289,6 +298,8 @@ using LinearAlgebra KernelGradient[i] = kernel_grad_acc ∇Cᵢ[i] = shift_c_acc ∇◌rᵢ[i] = shift_r_acc + ∇ρᵢ[i] = density_grad_acc + ∇uᵢ[i] = velocity_grad_acc AccelerationMax[i] = norm(acc_acc) end @@ -467,7 +478,7 @@ using LinearAlgebra SimMetaData::SimulationMetaData{D,T,S,K,B,L}, SimConstants, SimParticles, Position, Density, Pressure, Velocity, ParticleType, dρdt_acc, acc_acc, kernel_acc, kernel_grad_acc, shift_c_acc, - shift_r_acc, i, j) where {D,T,S<:ShiftingMode, + shift_r_acc, density_grad_acc, velocity_grad_acc, i, j) where {D,T,S<:ShiftingMode, K<:KernelOutputMode, B<:MDBCMode, L<:LogMode, @@ -512,16 +523,19 @@ using LinearAlgebra MotionLimiterCondition = MotionLimiterValue(eltype(ρᵢ), ParticleType[i]) * MotionLimiterValue(eltype(ρᵢ), ParticleType[j]) shift_c_acc += (m₀ / ρᵢ) * ∇ᵢWᵢⱼ shift_r_acc += (m₀ / ρⱼ) * dot(-xᵢⱼ, ∇ᵢWᵢⱼ) * MotionLimiterCondition + density_grad_acc += (m₀ / ρⱼ) * (ρⱼ - ρᵢ) * ∇ᵢWᵢⱼ + velocity_grad_acc += (m₀ / ρⱼ) * (vⱼ - vᵢ) * transpose(∇ᵢWᵢⱼ) end - return dρdt_acc, acc_acc, kernel_acc, kernel_grad_acc, shift_c_acc, shift_r_acc + return dρdt_acc, acc_acc, kernel_acc, kernel_grad_acc, shift_c_acc, shift_r_acc, + density_grad_acc, velocity_grad_acc end Base.@propagate_inbounds function ComputeInteractionsPerParticleNoKernel!( SimDensityDiffusion::SDD, SimViscosity::SV, SimKernel, SimMetaData::SimulationMetaData{D,T,S,NoKernelOutput,B,L}, SimConstants, SimParticles, Position, Density, Pressure, Velocity, ParticleType, - dρdt_acc, acc_acc, shift_c_acc, shift_r_acc, i, j) where {D,T, + dρdt_acc, acc_acc, shift_c_acc, shift_r_acc, density_grad_acc, velocity_grad_acc, i, j) where {D,T, S<:ShiftingMode, B<:MDBCMode, L<:LogMode, @@ -564,9 +578,11 @@ using LinearAlgebra MotionLimiterCondition = MotionLimiterValue(eltype(ρᵢ), ParticleType[i]) * MotionLimiterValue(eltype(ρᵢ), ParticleType[j]) shift_c_acc += (m₀ / ρᵢ) * ∇ᵢWᵢⱼ shift_r_acc += (m₀ / ρⱼ) * dot(-xᵢⱼ, ∇ᵢWᵢⱼ) * MotionLimiterCondition + density_grad_acc += (m₀ / ρⱼ) * (ρⱼ - ρᵢ) * ∇ᵢWᵢⱼ + velocity_grad_acc += (m₀ / ρⱼ) * (vⱼ - vᵢ) * transpose(∇ᵢWᵢⱼ) end - return dρdt_acc, acc_acc, shift_c_acc, shift_r_acc + return dρdt_acc, acc_acc, shift_c_acc, shift_r_acc, density_grad_acc, velocity_grad_acc end Base.@propagate_inbounds function ComputeInteractionsMDBC!(SimKernel, SimMetaData::SimulationMetaData{Dimensions, FloatType, SMode, KMode, BMode, LMode}, SimConstants, Position, Density, ParticleType, GhostPoints, i, j) where {Dimensions, FloatType, SMode, KMode, BMode, LMode} @@ -729,7 +745,7 @@ using LinearAlgebra @timeit SimMetaData.HourGlass "00b Init NeighborLoop" NeighborLoopPerParticle!( SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData, SimConstants, SimParticles, ParticleRanges, CellListIndices, - NeighborCellLists, dρdtI, SimParticles.Acceleration, ∇Cᵢ, ∇◌rᵢ, AccelerationMax, + NeighborCellLists, dρdtI, SimParticles.Acceleration, ∇Cᵢ, ∇◌rᵢ, ∇ρᵢ, ∇uᵢ, AccelerationMax, ) end @@ -761,11 +777,12 @@ using LinearAlgebra @timeit SimMetaData.HourGlass "02 Pressure" Pressure!(SimParticles.Pressure, SimParticles.Density, SimConstants) @timeit SimMetaData.HourGlass "03 Apply MDBC before Half TimeStep" ApplyMDBCBeforeHalf!(SimMetaData, SimKernel, SimConstants, SimParticles, ParticleRanges, UniqueCells) - @timeit SimMetaData.HourGlass "04 First NeighborLoop" NeighborLoopPerParticle!( - SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData, - SimConstants, SimParticles, ParticleRanges, CellListIndices, - NeighborCellLists, dρdtI, SimParticles.Acceleration, ∇Cᵢ, ∇◌rᵢ, AccelerationMax, - ) + @timeit SimMetaData.HourGlass "04 First NeighborLoop" NeighborLoopPerParticle!( + SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData, + SimConstants, SimParticles, ParticleRanges, CellListIndices, + NeighborCellLists, dρdtI, SimParticles.Acceleration, ∇Cᵢ, ∇◌rᵢ, + ∇ρᵢ, ∇uᵢ, AccelerationMax, + ) @timeit SimMetaData.HourGlass "05 Update To Half TimeStep" HalfTimeStep(SimMetaData, SimConstants, SimParticles, Positionₙ⁺, Velocityₙ⁺, ρₙ⁺, dρdtI, dt₂) @@ -777,7 +794,8 @@ using LinearAlgebra @timeit SimMetaData.HourGlass "08 Second NeighborLoop" NeighborLoopPerParticle!( SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData, SimConstants, SimParticles, ParticleRanges, CellListIndices, - NeighborCellLists, dρdtI, SimParticles.Acceleration, ∇Cᵢ, ∇◌rᵢ, AccelerationMax, + NeighborCellLists, dρdtI, SimParticles.Acceleration, ∇Cᵢ, ∇◌rᵢ, + ∇ρᵢ, ∇uᵢ, AccelerationMax, Position = Positionₙ⁺, Density = ρₙ⁺, Velocity = Velocityₙ⁺, @@ -795,7 +813,8 @@ using LinearAlgebra @timeit SimMetaData.HourGlass "06 NeighborLoop" NeighborLoopPerParticle!( SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData, SimConstants, SimParticles, ParticleRanges, CellListIndices, - NeighborCellLists, dρdtI, SimParticles.Acceleration, ∇Cᵢ, ∇◌rᵢ, AccelerationMax, + NeighborCellLists, dρdtI, SimParticles.Acceleration, ∇Cᵢ, ∇◌rᵢ, + ∇ρᵢ, ∇uᵢ, AccelerationMax, Position = Positionₙ⁺, Density = ρₙ⁺, Velocity = Velocityₙ⁺, @@ -806,7 +825,7 @@ using LinearAlgebra @timeit SimMetaData.HourGlass "08 Final LimitDensityAtBoundary" LimitDensityAtBoundary!(SimParticles.Density, SimConstants.ρ₀, ParticleType) - @timeit SimMetaData.HourGlass "09 Update To Final TimeStep" FullTimeStep(SimMetaData, SimKernel, SimConstants, SimParticles, Velocityₙ⁺, ∇Cᵢ, ∇◌rᵢ, dt) + @timeit SimMetaData.HourGlass "09 Update To Final TimeStep" FullTimeStep(SimMetaData, SimKernel, SimConstants, SimParticles, Velocityₙ⁺, ∇Cᵢ, ∇◌rᵢ, ∇ρᵢ, ∇uᵢ, dt) @timeit SimMetaData.HourGlass "10 Update MetaData" UpdateMetaData!(SimMetaData, dt) @@ -834,7 +853,8 @@ using LinearAlgebra SimMetaData.TimeSteppingMode = SimTimeStepping - dρdtI, Velocityₙ⁺, Positionₙ⁺, ρₙ⁺, ∇Cᵢ, ∇◌rᵢ = AllocateSupportDataStructures(SimMetaData, SimParticles.Position) + dρdtI, Velocityₙ⁺, Positionₙ⁺, ρₙ⁺, ∇Cᵢ, ∇◌rᵢ, ∇ρᵢ, ∇uᵢ = + AllocateSupportDataStructures(SimMetaData, SimParticles.Position) LoadMDBCNormals!(SimMetaData, SimParticles, ParticleNormalsPath) diff --git a/src/SPHExample.jl b/src/SPHExample.jl index 0ee0c0ae..2e381558 100644 --- a/src/SPHExample.jl +++ b/src/SPHExample.jl @@ -49,6 +49,7 @@ module SPHExample using .SimulationMetaDataConfiguration export SimulationMetaData, ShiftingMode, NoShifting, PlanarShifting, + FreeSurfaceMode, InternalFlow, FreeSurfaceCorrection, KernelOutputMode, NoKernelOutput, StoreKernelOutput, MDBCMode, NoMDBC, SimpleMDBC, LogMode, NoLog, StoreLog, diff --git a/src/SimulationMetaDataConfiguration.jl b/src/SimulationMetaDataConfiguration.jl index b92381bf..35a424ae 100644 --- a/src/SimulationMetaDataConfiguration.jl +++ b/src/SimulationMetaDataConfiguration.jl @@ -5,6 +5,7 @@ using TimerOutputs export SimulationMetaData, UpdateMetaData!, ShiftingMode, NoShifting, PlanarShifting, + FreeSurfaceMode, InternalFlow, FreeSurfaceCorrection, KernelOutputMode, NoKernelOutput, StoreKernelOutput, MDBCMode, NoMDBC, SimpleMDBC, LogMode, NoLog, StoreLog, @@ -14,6 +15,27 @@ abstract type ShiftingMode end struct NoShifting <: ShiftingMode end struct PlanarShifting <: ShiftingMode end +""" + FreeSurfaceMode + +Controls whether shifting uses a free-surface correction. +""" +abstract type FreeSurfaceMode end + +""" + InternalFlow + +Disables free-surface correction for internal flows. +""" +struct InternalFlow <: FreeSurfaceMode end + +""" + FreeSurfaceCorrection + +Enables free-surface correction constants for 2D/3D flows. +""" +struct FreeSurfaceCorrection <: FreeSurfaceMode end + abstract type KernelOutputMode end struct NoKernelOutput <: KernelOutputMode end struct StoreKernelOutput <: KernelOutputMode end @@ -55,6 +77,7 @@ struct SingleNeighborTimeStepping <: TimeSteppingMode end ExportGridCellParticleCounts::Bool = false OpenLogFile::Bool = true Δx::FloatType = zero(FloatType) + FreeSurfaceMode::FreeSurfaceMode = FreeSurfaceCorrection() TimeSteppingMode::TimeSteppingMode = SingleNeighborTimeStepping() end SimulationMetaData{D,T,S,K,B}(; kwargs...) where {D,T,S<:ShiftingMode,K<:KernelOutputMode,B<:MDBCMode} = diff --git a/src/TimeStepping.jl b/src/TimeStepping.jl index 4c2c4b9b..13aae93d 100644 --- a/src/TimeStepping.jl +++ b/src/TimeStepping.jl @@ -10,6 +10,18 @@ using ..SimulationEquations using ..SimulationGeometry using ..SimulationMetaDataConfiguration +@inline function FreeSurfaceThreshold(::Val{2}, ::Type{T}) where {T} + return T(1.5) +end + +@inline function FreeSurfaceThreshold(::Val{3}, ::Type{T}) where {T} + return T(2.75) +end + +@inline function FreeSurfaceThreshold(::Val{D}, ::Type{T}) where {D,T} + error("Free-surface shifting threshold is only defined for 2D and 3D.") +end + """ Δt(max_acceleration, SimulationConstants, SPHKernel) @@ -135,30 +147,48 @@ function FullTimeStep(::SimulationMetaData{D,T,NoShifting,K,B,L}, SimKernel, return nothing end +function FullTimeStep(SimMetaData::SimulationMetaData{D,T,NoShifting,K,B,L}, SimKernel, + SimConstants, SimParticles, Velocityₙ⁺, ∇Cᵢ, ∇◌rᵢ, _∇ρᵢ, _∇uᵢ, dt) where {D,T, + K<:KernelOutputMode, + B<:MDBCMode, + L<:LogMode} + return FullTimeStep(SimMetaData, SimKernel, SimConstants, SimParticles, Velocityₙ⁺, ∇Cᵢ, ∇◌rᵢ, dt) +end + function FullTimeStep(::SimulationMetaData{D,T,S,K,B,L}, SimKernel, SimConstants, - SimParticles, Velocityₙ⁺, ∇Cᵢ, ∇◌rᵢ, dt) where {D,T,S<:ShiftingMode, + SimParticles, Velocityₙ⁺, ∇Cᵢ, ∇◌rᵢ, ∇ρᵢ, ∇uᵢ, dt) where {D,T,S<:ShiftingMode, K<:KernelOutputMode, B<:MDBCMode, L<:LogMode} @unpack Position, Velocity, Acceleration = SimParticles + Density = SimParticles.Density ParticleType = SimParticles.Type AccelerationScalarType = eltype(eltype(Acceleration)) - A = 2# Value between 1 to 6 advised - A_FST = 0; # zero for internal flows - A_FSM = length(first(Position)); #2d, 3d val different + A = T(2) # Value between 1 to 6 advised + UseFreeSurfaceCorrection = SimMetaData.FreeSurfaceMode isa FreeSurfaceCorrection + A_FST = FreeSurfaceThreshold(Val(D), T) + A_FSM = T(D) @inbounds @simd ivdep for i in eachindex(Position) MotionLimiterFactor = MotionLimiterValue(AccelerationScalarType, ParticleType[i]) GravityFactor = GravityFactorValue(AccelerationScalarType, ParticleType[i]) Acceleration[i] += ConstructGravitySVector(Acceleration[i], SimConstants.g * GravityFactor) Velocity[i] += Acceleration[i] * dt * MotionLimiterFactor - A_FSC = (∇◌rᵢ[i] - A_FST)/(A_FSM - A_FST) - if A_FSC < 0 - δxᵢ = zero(eltype(Position)) + shift_scale = -A * SimKernel.h * norm(Velocityₙ⁺[i]) * dt + if UseFreeSurfaceCorrection + divergence_delta = ∇◌rᵢ[i] - A_FST + if divergence_delta < 0 + A_FSC = divergence_delta / (A_FSM - A_FST) + δxᵢ = shift_scale * A_FSC * ∇Cᵢ[i] + else + δxᵢ = shift_scale * ∇Cᵢ[i] + end else - δxᵢ = -A_FSC * A * SimKernel.h * norm(Velocityₙ⁺[i]) * dt * ∇Cᵢ[i] + δxᵢ = shift_scale * ∇Cᵢ[i] end + Density[i] += dot(δxᵢ, ∇ρᵢ[i]) * MotionLimiterFactor + Velocity[i] += (∇uᵢ[i] * δxᵢ) * MotionLimiterFactor Position[i] += (Velocityₙ⁺[i] * dt + δxᵢ) * MotionLimiterFactor end return nothing diff --git a/test/runtests.jl b/test/runtests.jl index 5ef7af40..87f63c41 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -46,7 +46,7 @@ end Type=typ, GroupMarker=group, GhostPoints=gpoint, GhostNormals=gnorm, )) - dρdtI, vel_n, pos_n, ρ_n, ∇C, ∇r = + dρdtI, vel_n, pos_n, ρ_n, ∇C, ∇r, ∇ρ, ∇u = AllocateSupportDataStructures(meta, particles.Position) for _ in 1:1000 From d420dc51fb6c34ae35a90b85d4f58161d7d8b46f Mon Sep 17 00:00:00 2001 From: AhmedSalih3d <36305327+AhmedSalih3d@users.noreply.github.com> Date: Mon, 26 Jan 2026 01:48:18 +0100 Subject: [PATCH 2/4] Add WCSPH shifting formula --- src/SPHCellList.jl | 117 +++++++++++++++++++++++++ src/SPHExample.jl | 2 +- src/SimulationMetaDataConfiguration.jl | 3 +- src/TimeStepping.jl | 30 +++++++ 4 files changed, 150 insertions(+), 2 deletions(-) diff --git a/src/SPHCellList.jl b/src/SPHCellList.jl index ff362d6b..e2fe8026 100644 --- a/src/SPHCellList.jl +++ b/src/SPHCellList.jl @@ -531,6 +531,67 @@ using LinearAlgebra density_grad_acc, velocity_grad_acc end + Base.@propagate_inbounds function ComputeInteractionsPerParticle!( + SimDensityDiffusion::SDD, SimViscosity::SV, SimKernel, + SimMetaData::SimulationMetaData{D,T,WCSPHShifting,K,B,L}, SimConstants, + SimParticles, Position, Density, Pressure, Velocity, ParticleType, + dρdt_acc, acc_acc, kernel_acc, kernel_grad_acc, shift_c_acc, + shift_r_acc, density_grad_acc, velocity_grad_acc, i, j) where {D,T, + K<:KernelOutputMode, + B<:MDBCMode, + L<:LogMode, + SDD<:SPHDensityDiffusion, + SV<:SPHViscosity} + @unpack m₀, dx = SimConstants + @unpack h⁻¹, H², h = SimKernel + + xᵢⱼ = Position[i] - Position[j] + xᵢⱼ² = dot(xᵢⱼ, xᵢⱼ) + if xᵢⱼ² <= H² + dᵢⱼ = sqrt(abs(xᵢⱼ²)) + dᵢⱼ² = dᵢⱼ^2 + q = clamp(dᵢⱼ * h⁻¹, 0.0, 2.0) + ∇ᵢWᵢⱼ = @fastpow ∇Wᵢⱼ(SimKernel, q, xᵢⱼ) + kernel_value = @fastpow Wᵢⱼ(SimKernel, q) + + ρᵢ = Density[i] + ρⱼ = Density[j] + + vᵢ = Velocity[i] + vⱼ = Velocity[j] + vᵢⱼ = vᵢ - vⱼ + density_symmetric_term = dot(-vᵢⱼ, ∇ᵢWᵢⱼ) + dρdt⁺ = -ρᵢ * (m₀ / ρⱼ) * density_symmetric_term + + Dᵢ, _ = compute_density_diffusion(SimDensityDiffusion, SimKernel, SimConstants, SimParticles, xᵢⱼ, ∇ᵢWᵢⱼ, dᵢⱼ², i, j, ParticleType) + + dρdt_acc += dρdt⁺ + Dᵢ + + Pᵢ = Pressure[i] + Pⱼ = Pressure[j] + Pfac = (Pᵢ + Pⱼ) / (ρᵢ * ρⱼ) + f_ab = tensile_correction(SimKernel, Pᵢ, ρᵢ, Pⱼ, ρⱼ, q, dx) + dvdt⁺ = -m₀ * (Pfac + f_ab) * ∇ᵢWᵢⱼ + + visc_term, _ = compute_viscosity(SimViscosity, SimKernel, SimConstants, SimParticles, xᵢⱼ, vᵢⱼ, ∇ᵢWᵢⱼ, dᵢⱼ², i, j) + + acc_acc += dvdt⁺ + visc_term + + kernel_acc, kernel_grad_acc = compute_kernel_output_local(SimMetaData, kernel_acc, kernel_grad_acc, SimKernel, q, ∇ᵢWᵢⱼ) + + w_ref = @fastpow Wᵢⱼ(SimKernel, dx * h⁻¹) + correction = one(T) + T(0.2) * (kernel_value / w_ref)^T(4) + shift_weight = (m₀ / (ρᵢ + ρⱼ)) * correction * kernel_value + shift_c_acc += shift_weight * xᵢⱼ + + density_grad_acc += (m₀ / ρⱼ) * (ρⱼ - ρᵢ) * ∇ᵢWᵢⱼ + velocity_grad_acc += (m₀ / ρⱼ) * (vⱼ - vᵢ) * transpose(∇ᵢWᵢⱼ) + end + + return dρdt_acc, acc_acc, kernel_acc, kernel_grad_acc, shift_c_acc, shift_r_acc, + density_grad_acc, velocity_grad_acc + end + Base.@propagate_inbounds function ComputeInteractionsPerParticleNoKernel!( SimDensityDiffusion::SDD, SimViscosity::SV, SimKernel, SimMetaData::SimulationMetaData{D,T,S,NoKernelOutput,B,L}, SimConstants, @@ -585,6 +646,62 @@ using LinearAlgebra return dρdt_acc, acc_acc, shift_c_acc, shift_r_acc, density_grad_acc, velocity_grad_acc end + Base.@propagate_inbounds function ComputeInteractionsPerParticleNoKernel!( + SimDensityDiffusion::SDD, SimViscosity::SV, SimKernel, + SimMetaData::SimulationMetaData{D,T,WCSPHShifting,NoKernelOutput,B,L}, SimConstants, + SimParticles, Position, Density, Pressure, Velocity, ParticleType, + dρdt_acc, acc_acc, shift_c_acc, shift_r_acc, density_grad_acc, velocity_grad_acc, i, j) where {D,T, + B<:MDBCMode, + L<:LogMode, + SDD<:SPHDensityDiffusion, + SV<:SPHViscosity} + @unpack m₀, dx = SimConstants + @unpack h⁻¹, H², h = SimKernel + + xᵢⱼ = Position[i] - Position[j] + xᵢⱼ² = dot(xᵢⱼ, xᵢⱼ) + if xᵢⱼ² <= H² + dᵢⱼ = sqrt(abs(xᵢⱼ²)) + dᵢⱼ² = dᵢⱼ^2 + q = clamp(dᵢⱼ * h⁻¹, 0.0, 2.0) + ∇ᵢWᵢⱼ = @fastpow ∇Wᵢⱼ(SimKernel, q, xᵢⱼ) + kernel_value = @fastpow Wᵢⱼ(SimKernel, q) + + ρᵢ = Density[i] + ρⱼ = Density[j] + + vᵢ = Velocity[i] + vⱼ = Velocity[j] + vᵢⱼ = vᵢ - vⱼ + density_symmetric_term = dot(-vᵢⱼ, ∇ᵢWᵢⱼ) + dρdt⁺ = -ρᵢ * (m₀ / ρⱼ) * density_symmetric_term + + Dᵢ, _ = compute_density_diffusion(SimDensityDiffusion, SimKernel, SimConstants, SimParticles, xᵢⱼ, ∇ᵢWᵢⱼ, dᵢⱼ², i, j, ParticleType) + + dρdt_acc += dρdt⁺ + Dᵢ + + Pᵢ = Pressure[i] + Pⱼ = Pressure[j] + Pfac = (Pᵢ + Pⱼ) / (ρᵢ * ρⱼ) + f_ab = tensile_correction(SimKernel, Pᵢ, ρᵢ, Pⱼ, ρⱼ, q, dx) + dvdt⁺ = -m₀ * (Pfac + f_ab) * ∇ᵢWᵢⱼ + + visc_term, _ = compute_viscosity(SimViscosity, SimKernel, SimConstants, SimParticles, xᵢⱼ, vᵢⱼ, ∇ᵢWᵢⱼ, dᵢⱼ², i, j) + + acc_acc += dvdt⁺ + visc_term + + w_ref = @fastpow Wᵢⱼ(SimKernel, dx * h⁻¹) + correction = one(T) + T(0.2) * (kernel_value / w_ref)^T(4) + shift_weight = (m₀ / (ρᵢ + ρⱼ)) * correction * kernel_value + shift_c_acc += shift_weight * xᵢⱼ + + density_grad_acc += (m₀ / ρⱼ) * (ρⱼ - ρᵢ) * ∇ᵢWᵢⱼ + velocity_grad_acc += (m₀ / ρⱼ) * (vⱼ - vᵢ) * transpose(∇ᵢWᵢⱼ) + end + + return dρdt_acc, acc_acc, shift_c_acc, shift_r_acc, density_grad_acc, velocity_grad_acc + end + Base.@propagate_inbounds function ComputeInteractionsMDBC!(SimKernel, SimMetaData::SimulationMetaData{Dimensions, FloatType, SMode, KMode, BMode, LMode}, SimConstants, Position, Density, ParticleType, GhostPoints, i, j) where {Dimensions, FloatType, SMode, KMode, BMode, LMode} @unpack ρ₀, m₀, α, γ, g, c₀, δᵩ, Cb, Cb⁻¹, ν₀, dx, SmagorinskyConstant, BlinConstant = SimConstants diff --git a/src/SPHExample.jl b/src/SPHExample.jl index 2e381558..9a95cf48 100644 --- a/src/SPHExample.jl +++ b/src/SPHExample.jl @@ -48,7 +48,7 @@ module SPHExample export SimulationLogger, generate_format_string, InitializeLogger, LogSimulationDetails, LogStep, LogFinal using .SimulationMetaDataConfiguration - export SimulationMetaData, ShiftingMode, NoShifting, PlanarShifting, + export SimulationMetaData, ShiftingMode, NoShifting, PlanarShifting, WCSPHShifting, FreeSurfaceMode, InternalFlow, FreeSurfaceCorrection, KernelOutputMode, NoKernelOutput, StoreKernelOutput, MDBCMode, NoMDBC, SimpleMDBC, diff --git a/src/SimulationMetaDataConfiguration.jl b/src/SimulationMetaDataConfiguration.jl index 35a424ae..0d3369cc 100644 --- a/src/SimulationMetaDataConfiguration.jl +++ b/src/SimulationMetaDataConfiguration.jl @@ -4,7 +4,7 @@ using Parameters using TimerOutputs -export SimulationMetaData, UpdateMetaData!, ShiftingMode, NoShifting, PlanarShifting, +export SimulationMetaData, UpdateMetaData!, ShiftingMode, NoShifting, PlanarShifting, WCSPHShifting, FreeSurfaceMode, InternalFlow, FreeSurfaceCorrection, KernelOutputMode, NoKernelOutput, StoreKernelOutput, MDBCMode, NoMDBC, SimpleMDBC, @@ -14,6 +14,7 @@ export SimulationMetaData, UpdateMetaData!, ShiftingMode, NoShifting, PlanarShif abstract type ShiftingMode end struct NoShifting <: ShiftingMode end struct PlanarShifting <: ShiftingMode end +struct WCSPHShifting <: ShiftingMode end """ FreeSurfaceMode diff --git a/src/TimeStepping.jl b/src/TimeStepping.jl index 13aae93d..60ddd310 100644 --- a/src/TimeStepping.jl +++ b/src/TimeStepping.jl @@ -194,4 +194,34 @@ function FullTimeStep(::SimulationMetaData{D,T,S,K,B,L}, SimKernel, SimConstants return nothing end +function FullTimeStep(::SimulationMetaData{D,T,WCSPHShifting,K,B,L}, SimKernel, SimConstants, + SimParticles, Velocityₙ⁺, ∇Cᵢ, _∇◌rᵢ, ∇ρᵢ, ∇uᵢ, dt) where {D,T, + K<:KernelOutputMode, + B<:MDBCMode, + L<:LogMode} + @unpack Position, Velocity, Acceleration = SimParticles + Density = SimParticles.Density + ParticleType = SimParticles.Type + AccelerationScalarType = eltype(eltype(Acceleration)) + + max_velocity = zero(eltype(eltype(Velocityₙ⁺))) + @inbounds for i in eachindex(Velocityₙ⁺) + max_velocity = max(max_velocity, norm(Velocityₙ⁺[i])) + end + shift_scale = -SimConstants.CFL * (max_velocity / SimConstants.c₀) * (2 * SimKernel.h)^2 + + @inbounds @simd ivdep for i in eachindex(Position) + MotionLimiterFactor = MotionLimiterValue(AccelerationScalarType, ParticleType[i]) + GravityFactor = GravityFactorValue(AccelerationScalarType, ParticleType[i]) + Acceleration[i] += ConstructGravitySVector(Acceleration[i], SimConstants.g * GravityFactor) + Velocity[i] += Acceleration[i] * dt * MotionLimiterFactor + + δxᵢ = shift_scale * ∇Cᵢ[i] + Density[i] += dot(δxᵢ, ∇ρᵢ[i]) * MotionLimiterFactor + Velocity[i] += (∇uᵢ[i] * δxᵢ) * MotionLimiterFactor + Position[i] += (Velocityₙ⁺[i] * dt + δxᵢ) * MotionLimiterFactor + end + return nothing +end + end From 24fd9e82541dc2227fad916dba73551dcc77689d Mon Sep 17 00:00:00 2001 From: AhmedSalih3d <36305327+AhmedSalih3d@users.noreply.github.com> Date: Mon, 26 Jan 2026 01:59:20 +0100 Subject: [PATCH 3/4] Implement WCSPH shifting equation --- src/PreProcess.jl | 8 +- src/SPHCellList.jl | 199 +++++-------------------- src/SPHExample.jl | 3 +- src/SimulationMetaDataConfiguration.jl | 26 +--- src/TimeStepping.jl | 69 +-------- test/runtests.jl | 2 +- 6 files changed, 42 insertions(+), 265 deletions(-) diff --git a/src/PreProcess.jl b/src/PreProcess.jl index ae98da7d..d00bc336 100644 --- a/src/PreProcess.jl +++ b/src/PreProcess.jl @@ -178,10 +178,8 @@ function AllocateSupportDataStructures(::SimulationMetaData{D,T,NoShifting,K,B,L ∇Cᵢ = Vector{PositionType}(undef, 0) ∇◌rᵢ = Vector{PositionUnderlyingType}(undef, 0) - ∇ρᵢ = Vector{PositionType}(undef, 0) - ∇uᵢ = Vector{SMatrix{D,D,PositionUnderlyingType}}(undef, 0) - return dρdtI, Velocityₙ⁺, Positionₙ⁺, ρₙ⁺, ∇Cᵢ, ∇◌rᵢ, ∇ρᵢ, ∇uᵢ + return dρdtI, Velocityₙ⁺, Positionₙ⁺, ρₙ⁺, ∇Cᵢ, ∇◌rᵢ end function AllocateSupportDataStructures(::SimulationMetaData{D,T,S,K,B,L}, Position) where {D,T,S<:ShiftingMode, @@ -200,10 +198,8 @@ function AllocateSupportDataStructures(::SimulationMetaData{D,T,S,K,B,L}, Positi ∇Cᵢ = zeros(PositionType, NumberOfPoints) ∇◌rᵢ = zeros(PositionUnderlyingType, NumberOfPoints) - ∇ρᵢ = zeros(PositionType, NumberOfPoints) - ∇uᵢ = zeros(SMatrix{D,D,PositionUnderlyingType}, NumberOfPoints) - return dρdtI, Velocityₙ⁺, Positionₙ⁺, ρₙ⁺, ∇Cᵢ, ∇◌rᵢ, ∇ρᵢ, ∇uᵢ + return dρdtI, Velocityₙ⁺, Positionₙ⁺, ρₙ⁺, ∇Cᵢ, ∇◌rᵢ end function LoadBoundaryNormals(::Val{D}, ::Type{T}, path_mdbc) where {D, T} diff --git a/src/SPHCellList.jl b/src/SPHCellList.jl index e2fe8026..2faddccb 100644 --- a/src/SPHCellList.jl +++ b/src/SPHCellList.jl @@ -37,7 +37,7 @@ using LinearAlgebra SimConstants, SimParticles, ParticleRanges, CellListIndices, NeighborCellLists, dρdtI, Acceleration, ∇Cᵢ, - ∇◌rᵢ, ∇ρᵢ, ∇uᵢ, AccelerationMax; + ∇◌rᵢ, AccelerationMax; Position = SimParticles.Position, Density = SimParticles.Density, Pressure = SimParticles.Pressure, @@ -93,7 +93,7 @@ using LinearAlgebra SimConstants, SimParticles, ParticleRanges, CellListIndices, NeighborCellLists, dρdtI, Acceleration, ∇Cᵢ, - ∇◌rᵢ, ∇ρᵢ, ∇uᵢ, AccelerationMax; + ∇◌rᵢ, AccelerationMax; Position = SimParticles.Position, Density = SimParticles.Density, Pressure = SimParticles.Pressure, @@ -161,7 +161,7 @@ using LinearAlgebra SimConstants, SimParticles, ParticleRanges, CellListIndices, NeighborCellLists, dρdtI, Acceleration, ∇Cᵢ, - ∇◌rᵢ, ∇ρᵢ, ∇uᵢ, AccelerationMax; + ∇◌rᵢ, AccelerationMax; Position = SimParticles.Position, Density = SimParticles.Density, Pressure = SimParticles.Pressure, @@ -175,41 +175,39 @@ using LinearAlgebra acc_acc = zero(Acceleration[i]) shift_c_acc = zero(∇Cᵢ[i]) shift_r_acc = zero(∇◌rᵢ[i]) - density_grad_acc = zero(∇ρᵢ[i]) - velocity_grad_acc = zero(∇uᵢ[i]) CellListIndex = CellListIndices[i] SameCellStart = ParticleRanges[CellListIndex] SameCellEnd = ParticleRanges[CellListIndex + 1] - 1 NeighborCellIndices = NeighborCellLists[CellListIndex] @inbounds for j in SameCellStart:(i - 1) - dρdt_acc, acc_acc, shift_c_acc, shift_r_acc, density_grad_acc, velocity_grad_acc = + dρdt_acc, acc_acc, shift_c_acc, shift_r_acc = ComputeInteractionsPerParticleNoKernel!( SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData, SimConstants, SimParticles, Position, Density, Pressure, Velocity, ParticleType, dρdt_acc, acc_acc, shift_c_acc, - shift_r_acc, density_grad_acc, velocity_grad_acc, i, j, + shift_r_acc, i, j, ) end @inbounds for j in (i + 1):SameCellEnd - dρdt_acc, acc_acc, shift_c_acc, shift_r_acc, density_grad_acc, velocity_grad_acc = + dρdt_acc, acc_acc, shift_c_acc, shift_r_acc = ComputeInteractionsPerParticleNoKernel!( SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData, SimConstants, SimParticles, Position, Density, Pressure, Velocity, ParticleType, dρdt_acc, acc_acc, shift_c_acc, - shift_r_acc, density_grad_acc, velocity_grad_acc, i, j, + shift_r_acc, i, j, ) end for NeighborIdx in NeighborCellIndices StartIndex_ = ParticleRanges[NeighborIdx] EndIndex_ = ParticleRanges[NeighborIdx + 1] - 1 @inbounds for j in StartIndex_:EndIndex_ - dρdt_acc, acc_acc, shift_c_acc, shift_r_acc, density_grad_acc, velocity_grad_acc = + dρdt_acc, acc_acc, shift_c_acc, shift_r_acc = ComputeInteractionsPerParticleNoKernel!( SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData, SimConstants, SimParticles, Position, Density, Pressure, Velocity, ParticleType, dρdt_acc, acc_acc, shift_c_acc, - shift_r_acc, density_grad_acc, velocity_grad_acc, i, j, + shift_r_acc, i, j, ) end end @@ -218,8 +216,6 @@ using LinearAlgebra Acceleration[i] = acc_acc ∇Cᵢ[i] = shift_c_acc ∇◌rᵢ[i] = shift_r_acc - ∇ρᵢ[i] = density_grad_acc - ∇uᵢ[i] = velocity_grad_acc AccelerationMax[i] = norm(acc_acc) end @@ -231,7 +227,7 @@ using LinearAlgebra SimConstants, SimParticles, ParticleRanges, CellListIndices, NeighborCellLists, dρdtI, Acceleration, ∇Cᵢ, - ∇◌rᵢ, ∇ρᵢ, ∇uᵢ, AccelerationMax; + ∇◌rᵢ, AccelerationMax; Position = SimParticles.Position, Density = SimParticles.Density, Pressure = SimParticles.Pressure, @@ -250,8 +246,6 @@ using LinearAlgebra kernel_grad_acc = zero(KernelGradient[i]) shift_c_acc = zero(∇Cᵢ[i]) shift_r_acc = zero(∇◌rᵢ[i]) - density_grad_acc = zero(∇ρᵢ[i]) - velocity_grad_acc = zero(∇uᵢ[i]) CellListIndex = CellListIndices[i] SameCellStart = ParticleRanges[CellListIndex] SameCellEnd = ParticleRanges[CellListIndex + 1] - 1 @@ -259,22 +253,20 @@ using LinearAlgebra @inbounds for j in SameCellStart:(i - 1) dρdt_acc, acc_acc, kernel_acc, kernel_grad_acc, shift_c_acc, - shift_r_acc, density_grad_acc, velocity_grad_acc = ComputeInteractionsPerParticle!( + shift_r_acc = ComputeInteractionsPerParticle!( SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData, SimConstants, SimParticles, Position, Density, Pressure, Velocity, ParticleType, dρdt_acc, acc_acc, kernel_acc, - kernel_grad_acc, shift_c_acc, shift_r_acc, density_grad_acc, - velocity_grad_acc, i, j, + kernel_grad_acc, shift_c_acc, shift_r_acc, i, j, ) end @inbounds for j in (i + 1):SameCellEnd dρdt_acc, acc_acc, kernel_acc, kernel_grad_acc, shift_c_acc, - shift_r_acc, density_grad_acc, velocity_grad_acc = ComputeInteractionsPerParticle!( + shift_r_acc = ComputeInteractionsPerParticle!( SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData, SimConstants, SimParticles, Position, Density, Pressure, Velocity, ParticleType, dρdt_acc, acc_acc, kernel_acc, - kernel_grad_acc, shift_c_acc, shift_r_acc, density_grad_acc, - velocity_grad_acc, i, j, + kernel_grad_acc, shift_c_acc, shift_r_acc, i, j, ) end for NeighborIdx in NeighborCellIndices @@ -282,12 +274,11 @@ using LinearAlgebra EndIndex_ = ParticleRanges[NeighborIdx + 1] - 1 @inbounds for j in StartIndex_:EndIndex_ dρdt_acc, acc_acc, kernel_acc, kernel_grad_acc, shift_c_acc, - shift_r_acc, density_grad_acc, velocity_grad_acc = ComputeInteractionsPerParticle!( + shift_r_acc = ComputeInteractionsPerParticle!( SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData, SimConstants, SimParticles, Position, Density, Pressure, Velocity, ParticleType, dρdt_acc, acc_acc, kernel_acc, - kernel_grad_acc, shift_c_acc, shift_r_acc, density_grad_acc, - velocity_grad_acc, i, j, + kernel_grad_acc, shift_c_acc, shift_r_acc, i, j, ) end end @@ -298,8 +289,6 @@ using LinearAlgebra KernelGradient[i] = kernel_grad_acc ∇Cᵢ[i] = shift_c_acc ∇◌rᵢ[i] = shift_r_acc - ∇ρᵢ[i] = density_grad_acc - ∇uᵢ[i] = velocity_grad_acc AccelerationMax[i] = norm(acc_acc) end @@ -478,7 +467,7 @@ using LinearAlgebra SimMetaData::SimulationMetaData{D,T,S,K,B,L}, SimConstants, SimParticles, Position, Density, Pressure, Velocity, ParticleType, dρdt_acc, acc_acc, kernel_acc, kernel_grad_acc, shift_c_acc, - shift_r_acc, density_grad_acc, velocity_grad_acc, i, j) where {D,T,S<:ShiftingMode, + shift_r_acc, i, j) where {D,T,S<:ShiftingMode, K<:KernelOutputMode, B<:MDBCMode, L<:LogMode, @@ -521,82 +510,21 @@ using LinearAlgebra kernel_acc, kernel_grad_acc = compute_kernel_output_local(SimMetaData, kernel_acc, kernel_grad_acc, SimKernel, q, ∇ᵢWᵢⱼ) MotionLimiterCondition = MotionLimiterValue(eltype(ρᵢ), ParticleType[i]) * MotionLimiterValue(eltype(ρᵢ), ParticleType[j]) - shift_c_acc += (m₀ / ρᵢ) * ∇ᵢWᵢⱼ - shift_r_acc += (m₀ / ρⱼ) * dot(-xᵢⱼ, ∇ᵢWᵢⱼ) * MotionLimiterCondition - density_grad_acc += (m₀ / ρⱼ) * (ρⱼ - ρᵢ) * ∇ᵢWᵢⱼ - velocity_grad_acc += (m₀ / ρⱼ) * (vⱼ - vᵢ) * transpose(∇ᵢWᵢⱼ) - end - - return dρdt_acc, acc_acc, kernel_acc, kernel_grad_acc, shift_c_acc, shift_r_acc, - density_grad_acc, velocity_grad_acc - end - - Base.@propagate_inbounds function ComputeInteractionsPerParticle!( - SimDensityDiffusion::SDD, SimViscosity::SV, SimKernel, - SimMetaData::SimulationMetaData{D,T,WCSPHShifting,K,B,L}, SimConstants, - SimParticles, Position, Density, Pressure, Velocity, ParticleType, - dρdt_acc, acc_acc, kernel_acc, kernel_grad_acc, shift_c_acc, - shift_r_acc, density_grad_acc, velocity_grad_acc, i, j) where {D,T, - K<:KernelOutputMode, - B<:MDBCMode, - L<:LogMode, - SDD<:SPHDensityDiffusion, - SV<:SPHViscosity} - @unpack m₀, dx = SimConstants - @unpack h⁻¹, H², h = SimKernel - - xᵢⱼ = Position[i] - Position[j] - xᵢⱼ² = dot(xᵢⱼ, xᵢⱼ) - if xᵢⱼ² <= H² - dᵢⱼ = sqrt(abs(xᵢⱼ²)) - dᵢⱼ² = dᵢⱼ^2 - q = clamp(dᵢⱼ * h⁻¹, 0.0, 2.0) - ∇ᵢWᵢⱼ = @fastpow ∇Wᵢⱼ(SimKernel, q, xᵢⱼ) kernel_value = @fastpow Wᵢⱼ(SimKernel, q) - - ρᵢ = Density[i] - ρⱼ = Density[j] - - vᵢ = Velocity[i] - vⱼ = Velocity[j] - vᵢⱼ = vᵢ - vⱼ - density_symmetric_term = dot(-vᵢⱼ, ∇ᵢWᵢⱼ) - dρdt⁺ = -ρᵢ * (m₀ / ρⱼ) * density_symmetric_term - - Dᵢ, _ = compute_density_diffusion(SimDensityDiffusion, SimKernel, SimConstants, SimParticles, xᵢⱼ, ∇ᵢWᵢⱼ, dᵢⱼ², i, j, ParticleType) - - dρdt_acc += dρdt⁺ + Dᵢ - - Pᵢ = Pressure[i] - Pⱼ = Pressure[j] - Pfac = (Pᵢ + Pⱼ) / (ρᵢ * ρⱼ) - f_ab = tensile_correction(SimKernel, Pᵢ, ρᵢ, Pⱼ, ρⱼ, q, dx) - dvdt⁺ = -m₀ * (Pfac + f_ab) * ∇ᵢWᵢⱼ - - visc_term, _ = compute_viscosity(SimViscosity, SimKernel, SimConstants, SimParticles, xᵢⱼ, vᵢⱼ, ∇ᵢWᵢⱼ, dᵢⱼ², i, j) - - acc_acc += dvdt⁺ + visc_term - - kernel_acc, kernel_grad_acc = compute_kernel_output_local(SimMetaData, kernel_acc, kernel_grad_acc, SimKernel, q, ∇ᵢWᵢⱼ) - w_ref = @fastpow Wᵢⱼ(SimKernel, dx * h⁻¹) correction = one(T) + T(0.2) * (kernel_value / w_ref)^T(4) - shift_weight = (m₀ / (ρᵢ + ρⱼ)) * correction * kernel_value - shift_c_acc += shift_weight * xᵢⱼ - - density_grad_acc += (m₀ / ρⱼ) * (ρⱼ - ρᵢ) * ∇ᵢWᵢⱼ - velocity_grad_acc += (m₀ / ρⱼ) * (vⱼ - vᵢ) * transpose(∇ᵢWᵢⱼ) + shift_weight = (m₀ / ρᵢ) * (m₀ / (ρᵢ + ρⱼ)) * correction * kernel_value + shift_c_acc += shift_weight * xᵢⱼ * MotionLimiterCondition end - return dρdt_acc, acc_acc, kernel_acc, kernel_grad_acc, shift_c_acc, shift_r_acc, - density_grad_acc, velocity_grad_acc + return dρdt_acc, acc_acc, kernel_acc, kernel_grad_acc, shift_c_acc, shift_r_acc end Base.@propagate_inbounds function ComputeInteractionsPerParticleNoKernel!( SimDensityDiffusion::SDD, SimViscosity::SV, SimKernel, SimMetaData::SimulationMetaData{D,T,S,NoKernelOutput,B,L}, SimConstants, SimParticles, Position, Density, Pressure, Velocity, ParticleType, - dρdt_acc, acc_acc, shift_c_acc, shift_r_acc, density_grad_acc, velocity_grad_acc, i, j) where {D,T, + dρdt_acc, acc_acc, shift_c_acc, shift_r_acc, i, j) where {D,T, S<:ShiftingMode, B<:MDBCMode, L<:LogMode, @@ -637,69 +565,14 @@ using LinearAlgebra acc_acc += dvdt⁺ + visc_term MotionLimiterCondition = MotionLimiterValue(eltype(ρᵢ), ParticleType[i]) * MotionLimiterValue(eltype(ρᵢ), ParticleType[j]) - shift_c_acc += (m₀ / ρᵢ) * ∇ᵢWᵢⱼ - shift_r_acc += (m₀ / ρⱼ) * dot(-xᵢⱼ, ∇ᵢWᵢⱼ) * MotionLimiterCondition - density_grad_acc += (m₀ / ρⱼ) * (ρⱼ - ρᵢ) * ∇ᵢWᵢⱼ - velocity_grad_acc += (m₀ / ρⱼ) * (vⱼ - vᵢ) * transpose(∇ᵢWᵢⱼ) - end - - return dρdt_acc, acc_acc, shift_c_acc, shift_r_acc, density_grad_acc, velocity_grad_acc - end - - Base.@propagate_inbounds function ComputeInteractionsPerParticleNoKernel!( - SimDensityDiffusion::SDD, SimViscosity::SV, SimKernel, - SimMetaData::SimulationMetaData{D,T,WCSPHShifting,NoKernelOutput,B,L}, SimConstants, - SimParticles, Position, Density, Pressure, Velocity, ParticleType, - dρdt_acc, acc_acc, shift_c_acc, shift_r_acc, density_grad_acc, velocity_grad_acc, i, j) where {D,T, - B<:MDBCMode, - L<:LogMode, - SDD<:SPHDensityDiffusion, - SV<:SPHViscosity} - @unpack m₀, dx = SimConstants - @unpack h⁻¹, H², h = SimKernel - - xᵢⱼ = Position[i] - Position[j] - xᵢⱼ² = dot(xᵢⱼ, xᵢⱼ) - if xᵢⱼ² <= H² - dᵢⱼ = sqrt(abs(xᵢⱼ²)) - dᵢⱼ² = dᵢⱼ^2 - q = clamp(dᵢⱼ * h⁻¹, 0.0, 2.0) - ∇ᵢWᵢⱼ = @fastpow ∇Wᵢⱼ(SimKernel, q, xᵢⱼ) kernel_value = @fastpow Wᵢⱼ(SimKernel, q) - - ρᵢ = Density[i] - ρⱼ = Density[j] - - vᵢ = Velocity[i] - vⱼ = Velocity[j] - vᵢⱼ = vᵢ - vⱼ - density_symmetric_term = dot(-vᵢⱼ, ∇ᵢWᵢⱼ) - dρdt⁺ = -ρᵢ * (m₀ / ρⱼ) * density_symmetric_term - - Dᵢ, _ = compute_density_diffusion(SimDensityDiffusion, SimKernel, SimConstants, SimParticles, xᵢⱼ, ∇ᵢWᵢⱼ, dᵢⱼ², i, j, ParticleType) - - dρdt_acc += dρdt⁺ + Dᵢ - - Pᵢ = Pressure[i] - Pⱼ = Pressure[j] - Pfac = (Pᵢ + Pⱼ) / (ρᵢ * ρⱼ) - f_ab = tensile_correction(SimKernel, Pᵢ, ρᵢ, Pⱼ, ρⱼ, q, dx) - dvdt⁺ = -m₀ * (Pfac + f_ab) * ∇ᵢWᵢⱼ - - visc_term, _ = compute_viscosity(SimViscosity, SimKernel, SimConstants, SimParticles, xᵢⱼ, vᵢⱼ, ∇ᵢWᵢⱼ, dᵢⱼ², i, j) - - acc_acc += dvdt⁺ + visc_term - w_ref = @fastpow Wᵢⱼ(SimKernel, dx * h⁻¹) correction = one(T) + T(0.2) * (kernel_value / w_ref)^T(4) - shift_weight = (m₀ / (ρᵢ + ρⱼ)) * correction * kernel_value - shift_c_acc += shift_weight * xᵢⱼ - - density_grad_acc += (m₀ / ρⱼ) * (ρⱼ - ρᵢ) * ∇ᵢWᵢⱼ - velocity_grad_acc += (m₀ / ρⱼ) * (vⱼ - vᵢ) * transpose(∇ᵢWᵢⱼ) + shift_weight = (m₀ / ρᵢ) * (m₀ / (ρᵢ + ρⱼ)) * correction * kernel_value + shift_c_acc += shift_weight * xᵢⱼ * MotionLimiterCondition end - return dρdt_acc, acc_acc, shift_c_acc, shift_r_acc, density_grad_acc, velocity_grad_acc + return dρdt_acc, acc_acc, shift_c_acc, shift_r_acc end Base.@propagate_inbounds function ComputeInteractionsMDBC!(SimKernel, SimMetaData::SimulationMetaData{Dimensions, FloatType, SMode, KMode, BMode, LMode}, SimConstants, Position, Density, ParticleType, GhostPoints, i, j) where {Dimensions, FloatType, SMode, KMode, BMode, LMode} @@ -862,7 +735,7 @@ using LinearAlgebra @timeit SimMetaData.HourGlass "00b Init NeighborLoop" NeighborLoopPerParticle!( SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData, SimConstants, SimParticles, ParticleRanges, CellListIndices, - NeighborCellLists, dρdtI, SimParticles.Acceleration, ∇Cᵢ, ∇◌rᵢ, ∇ρᵢ, ∇uᵢ, AccelerationMax, + NeighborCellLists, dρdtI, SimParticles.Acceleration, ∇Cᵢ, ∇◌rᵢ, AccelerationMax, ) end @@ -894,12 +767,11 @@ using LinearAlgebra @timeit SimMetaData.HourGlass "02 Pressure" Pressure!(SimParticles.Pressure, SimParticles.Density, SimConstants) @timeit SimMetaData.HourGlass "03 Apply MDBC before Half TimeStep" ApplyMDBCBeforeHalf!(SimMetaData, SimKernel, SimConstants, SimParticles, ParticleRanges, UniqueCells) - @timeit SimMetaData.HourGlass "04 First NeighborLoop" NeighborLoopPerParticle!( - SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData, - SimConstants, SimParticles, ParticleRanges, CellListIndices, - NeighborCellLists, dρdtI, SimParticles.Acceleration, ∇Cᵢ, ∇◌rᵢ, - ∇ρᵢ, ∇uᵢ, AccelerationMax, - ) + @timeit SimMetaData.HourGlass "04 First NeighborLoop" NeighborLoopPerParticle!( + SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData, + SimConstants, SimParticles, ParticleRanges, CellListIndices, + NeighborCellLists, dρdtI, SimParticles.Acceleration, ∇Cᵢ, ∇◌rᵢ, AccelerationMax, + ) @timeit SimMetaData.HourGlass "05 Update To Half TimeStep" HalfTimeStep(SimMetaData, SimConstants, SimParticles, Positionₙ⁺, Velocityₙ⁺, ρₙ⁺, dρdtI, dt₂) @@ -911,8 +783,7 @@ using LinearAlgebra @timeit SimMetaData.HourGlass "08 Second NeighborLoop" NeighborLoopPerParticle!( SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData, SimConstants, SimParticles, ParticleRanges, CellListIndices, - NeighborCellLists, dρdtI, SimParticles.Acceleration, ∇Cᵢ, ∇◌rᵢ, - ∇ρᵢ, ∇uᵢ, AccelerationMax, + NeighborCellLists, dρdtI, SimParticles.Acceleration, ∇Cᵢ, ∇◌rᵢ, AccelerationMax, Position = Positionₙ⁺, Density = ρₙ⁺, Velocity = Velocityₙ⁺, @@ -930,8 +801,7 @@ using LinearAlgebra @timeit SimMetaData.HourGlass "06 NeighborLoop" NeighborLoopPerParticle!( SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData, SimConstants, SimParticles, ParticleRanges, CellListIndices, - NeighborCellLists, dρdtI, SimParticles.Acceleration, ∇Cᵢ, ∇◌rᵢ, - ∇ρᵢ, ∇uᵢ, AccelerationMax, + NeighborCellLists, dρdtI, SimParticles.Acceleration, ∇Cᵢ, ∇◌rᵢ, AccelerationMax, Position = Positionₙ⁺, Density = ρₙ⁺, Velocity = Velocityₙ⁺, @@ -942,7 +812,7 @@ using LinearAlgebra @timeit SimMetaData.HourGlass "08 Final LimitDensityAtBoundary" LimitDensityAtBoundary!(SimParticles.Density, SimConstants.ρ₀, ParticleType) - @timeit SimMetaData.HourGlass "09 Update To Final TimeStep" FullTimeStep(SimMetaData, SimKernel, SimConstants, SimParticles, Velocityₙ⁺, ∇Cᵢ, ∇◌rᵢ, ∇ρᵢ, ∇uᵢ, dt) + @timeit SimMetaData.HourGlass "09 Update To Final TimeStep" FullTimeStep(SimMetaData, SimKernel, SimConstants, SimParticles, Velocityₙ⁺, ∇Cᵢ, ∇◌rᵢ, dt) @timeit SimMetaData.HourGlass "10 Update MetaData" UpdateMetaData!(SimMetaData, dt) @@ -970,8 +840,7 @@ using LinearAlgebra SimMetaData.TimeSteppingMode = SimTimeStepping - dρdtI, Velocityₙ⁺, Positionₙ⁺, ρₙ⁺, ∇Cᵢ, ∇◌rᵢ, ∇ρᵢ, ∇uᵢ = - AllocateSupportDataStructures(SimMetaData, SimParticles.Position) + dρdtI, Velocityₙ⁺, Positionₙ⁺, ρₙ⁺, ∇Cᵢ, ∇◌rᵢ = AllocateSupportDataStructures(SimMetaData, SimParticles.Position) LoadMDBCNormals!(SimMetaData, SimParticles, ParticleNormalsPath) diff --git a/src/SPHExample.jl b/src/SPHExample.jl index 9a95cf48..0ee0c0ae 100644 --- a/src/SPHExample.jl +++ b/src/SPHExample.jl @@ -48,8 +48,7 @@ module SPHExample export SimulationLogger, generate_format_string, InitializeLogger, LogSimulationDetails, LogStep, LogFinal using .SimulationMetaDataConfiguration - export SimulationMetaData, ShiftingMode, NoShifting, PlanarShifting, WCSPHShifting, - FreeSurfaceMode, InternalFlow, FreeSurfaceCorrection, + export SimulationMetaData, ShiftingMode, NoShifting, PlanarShifting, KernelOutputMode, NoKernelOutput, StoreKernelOutput, MDBCMode, NoMDBC, SimpleMDBC, LogMode, NoLog, StoreLog, diff --git a/src/SimulationMetaDataConfiguration.jl b/src/SimulationMetaDataConfiguration.jl index 0d3369cc..b92381bf 100644 --- a/src/SimulationMetaDataConfiguration.jl +++ b/src/SimulationMetaDataConfiguration.jl @@ -4,8 +4,7 @@ using Parameters using TimerOutputs -export SimulationMetaData, UpdateMetaData!, ShiftingMode, NoShifting, PlanarShifting, WCSPHShifting, - FreeSurfaceMode, InternalFlow, FreeSurfaceCorrection, +export SimulationMetaData, UpdateMetaData!, ShiftingMode, NoShifting, PlanarShifting, KernelOutputMode, NoKernelOutput, StoreKernelOutput, MDBCMode, NoMDBC, SimpleMDBC, LogMode, NoLog, StoreLog, @@ -14,28 +13,6 @@ export SimulationMetaData, UpdateMetaData!, ShiftingMode, NoShifting, PlanarShif abstract type ShiftingMode end struct NoShifting <: ShiftingMode end struct PlanarShifting <: ShiftingMode end -struct WCSPHShifting <: ShiftingMode end - -""" - FreeSurfaceMode - -Controls whether shifting uses a free-surface correction. -""" -abstract type FreeSurfaceMode end - -""" - InternalFlow - -Disables free-surface correction for internal flows. -""" -struct InternalFlow <: FreeSurfaceMode end - -""" - FreeSurfaceCorrection - -Enables free-surface correction constants for 2D/3D flows. -""" -struct FreeSurfaceCorrection <: FreeSurfaceMode end abstract type KernelOutputMode end struct NoKernelOutput <: KernelOutputMode end @@ -78,7 +55,6 @@ struct SingleNeighborTimeStepping <: TimeSteppingMode end ExportGridCellParticleCounts::Bool = false OpenLogFile::Bool = true Δx::FloatType = zero(FloatType) - FreeSurfaceMode::FreeSurfaceMode = FreeSurfaceCorrection() TimeSteppingMode::TimeSteppingMode = SingleNeighborTimeStepping() end SimulationMetaData{D,T,S,K,B}(; kwargs...) where {D,T,S<:ShiftingMode,K<:KernelOutputMode,B<:MDBCMode} = diff --git a/src/TimeStepping.jl b/src/TimeStepping.jl index 60ddd310..cbbda82a 100644 --- a/src/TimeStepping.jl +++ b/src/TimeStepping.jl @@ -10,18 +10,6 @@ using ..SimulationEquations using ..SimulationGeometry using ..SimulationMetaDataConfiguration -@inline function FreeSurfaceThreshold(::Val{2}, ::Type{T}) where {T} - return T(1.5) -end - -@inline function FreeSurfaceThreshold(::Val{3}, ::Type{T}) where {T} - return T(2.75) -end - -@inline function FreeSurfaceThreshold(::Val{D}, ::Type{T}) where {D,T} - error("Free-surface shifting threshold is only defined for 2D and 3D.") -end - """ Δt(max_acceleration, SimulationConstants, SPHKernel) @@ -147,69 +135,19 @@ function FullTimeStep(::SimulationMetaData{D,T,NoShifting,K,B,L}, SimKernel, return nothing end -function FullTimeStep(SimMetaData::SimulationMetaData{D,T,NoShifting,K,B,L}, SimKernel, - SimConstants, SimParticles, Velocityₙ⁺, ∇Cᵢ, ∇◌rᵢ, _∇ρᵢ, _∇uᵢ, dt) where {D,T, - K<:KernelOutputMode, - B<:MDBCMode, - L<:LogMode} - return FullTimeStep(SimMetaData, SimKernel, SimConstants, SimParticles, Velocityₙ⁺, ∇Cᵢ, ∇◌rᵢ, dt) -end - function FullTimeStep(::SimulationMetaData{D,T,S,K,B,L}, SimKernel, SimConstants, - SimParticles, Velocityₙ⁺, ∇Cᵢ, ∇◌rᵢ, ∇ρᵢ, ∇uᵢ, dt) where {D,T,S<:ShiftingMode, + SimParticles, Velocityₙ⁺, ∇Cᵢ, ∇◌rᵢ, dt) where {D,T,S<:ShiftingMode, K<:KernelOutputMode, B<:MDBCMode, L<:LogMode} @unpack Position, Velocity, Acceleration = SimParticles - Density = SimParticles.Density ParticleType = SimParticles.Type AccelerationScalarType = eltype(eltype(Acceleration)) - A = T(2) # Value between 1 to 6 advised - UseFreeSurfaceCorrection = SimMetaData.FreeSurfaceMode isa FreeSurfaceCorrection - A_FST = FreeSurfaceThreshold(Val(D), T) - A_FSM = T(D) - @inbounds @simd ivdep for i in eachindex(Position) - MotionLimiterFactor = MotionLimiterValue(AccelerationScalarType, ParticleType[i]) - GravityFactor = GravityFactorValue(AccelerationScalarType, ParticleType[i]) - Acceleration[i] += ConstructGravitySVector(Acceleration[i], SimConstants.g * GravityFactor) - Velocity[i] += Acceleration[i] * dt * MotionLimiterFactor - - shift_scale = -A * SimKernel.h * norm(Velocityₙ⁺[i]) * dt - if UseFreeSurfaceCorrection - divergence_delta = ∇◌rᵢ[i] - A_FST - if divergence_delta < 0 - A_FSC = divergence_delta / (A_FSM - A_FST) - δxᵢ = shift_scale * A_FSC * ∇Cᵢ[i] - else - δxᵢ = shift_scale * ∇Cᵢ[i] - end - else - δxᵢ = shift_scale * ∇Cᵢ[i] - end - - Density[i] += dot(δxᵢ, ∇ρᵢ[i]) * MotionLimiterFactor - Velocity[i] += (∇uᵢ[i] * δxᵢ) * MotionLimiterFactor - Position[i] += (Velocityₙ⁺[i] * dt + δxᵢ) * MotionLimiterFactor - end - return nothing -end - -function FullTimeStep(::SimulationMetaData{D,T,WCSPHShifting,K,B,L}, SimKernel, SimConstants, - SimParticles, Velocityₙ⁺, ∇Cᵢ, _∇◌rᵢ, ∇ρᵢ, ∇uᵢ, dt) where {D,T, - K<:KernelOutputMode, - B<:MDBCMode, - L<:LogMode} - @unpack Position, Velocity, Acceleration = SimParticles - Density = SimParticles.Density - ParticleType = SimParticles.Type - AccelerationScalarType = eltype(eltype(Acceleration)) - max_velocity = zero(eltype(eltype(Velocityₙ⁺))) @inbounds for i in eachindex(Velocityₙ⁺) max_velocity = max(max_velocity, norm(Velocityₙ⁺[i])) end shift_scale = -SimConstants.CFL * (max_velocity / SimConstants.c₀) * (2 * SimKernel.h)^2 - @inbounds @simd ivdep for i in eachindex(Position) MotionLimiterFactor = MotionLimiterValue(AccelerationScalarType, ParticleType[i]) GravityFactor = GravityFactorValue(AccelerationScalarType, ParticleType[i]) @@ -217,9 +155,8 @@ function FullTimeStep(::SimulationMetaData{D,T,WCSPHShifting,K,B,L}, SimKernel, Velocity[i] += Acceleration[i] * dt * MotionLimiterFactor δxᵢ = shift_scale * ∇Cᵢ[i] - Density[i] += dot(δxᵢ, ∇ρᵢ[i]) * MotionLimiterFactor - Velocity[i] += (∇uᵢ[i] * δxᵢ) * MotionLimiterFactor - Position[i] += (Velocityₙ⁺[i] * dt + δxᵢ) * MotionLimiterFactor + + Position[i] += (Velocityₙ⁺[i] * dt + δxᵢ) * MotionLimiterFactor end return nothing end diff --git a/test/runtests.jl b/test/runtests.jl index 87f63c41..5ef7af40 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -46,7 +46,7 @@ end Type=typ, GroupMarker=group, GhostPoints=gpoint, GhostNormals=gnorm, )) - dρdtI, vel_n, pos_n, ρ_n, ∇C, ∇r, ∇ρ, ∇u = + dρdtI, vel_n, pos_n, ρ_n, ∇C, ∇r = AllocateSupportDataStructures(meta, particles.Position) for _ in 1:1000 From 6c3b8ff42df469107db58ce30c87475167db63cd Mon Sep 17 00:00:00 2001 From: AhmedSalih3d <36305327+AhmedSalih3d@users.noreply.github.com> Date: Mon, 26 Jan 2026 02:32:42 +0100 Subject: [PATCH 4/4] Add shift-specific CFL --- src/SimulationConstantsConfiguration.jl | 2 ++ src/TimeStepping.jl | 2 +- 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/src/SimulationConstantsConfiguration.jl b/src/SimulationConstantsConfiguration.jl index 7d3b569f..0d6a5fe6 100644 --- a/src/SimulationConstantsConfiguration.jl +++ b/src/SimulationConstantsConfiguration.jl @@ -23,6 +23,7 @@ SimulationConstants is a parameterized struct representing the constants and par - `dt::T`: Initial time step. Default is 1e-5. - `δᵩ::T`: Coefficient for density diffusion. Default is 0.1. - `CFL::T`: CFL (Courant-Friedrichs-Lewy) number (positive). Default is 0.2. +- `ShiftCFL::T`: CFL value used for shifting (positive). Default is 0.2. - `η²::T`: Eta squared (positive). Default is computed as `(0.01 * H)^2`. # Example @@ -44,6 +45,7 @@ constants = SimulationConstants(ρ₀=1017, dx=0.03, α=0.02) γ⁻¹::T = 1/γ ; @assert γ⁻¹ > 0 "Inverse adiabatic index (γ⁻¹) must be positive" δᵩ::T = 0.1 ; @assert δᵩ > 0 "Density variation (δᵩ) must be positive" CFL::T = 0.2 ; @assert CFL > 0 "CFL condition (CFL) must be positive" + ShiftCFL::T = 0.2 ; @assert ShiftCFL > 0 "Shift CFL condition (ShiftCFL) must be positive" Cb::T = (c₀^2 * ρ₀)/γ ; @assert Cb >= 0 "Cb (pressure coefficient) must be positive" Cb⁻¹::T = inv(Cb) ; @assert Cb⁻¹>= 0 "Inverse Cb (inverse pressure coefficient) must be positive" ν₀::T = 1e-6 ; @assert ν₀ >= 0 "Kinematic viscosity must be positive" diff --git a/src/TimeStepping.jl b/src/TimeStepping.jl index cbbda82a..84074a5e 100644 --- a/src/TimeStepping.jl +++ b/src/TimeStepping.jl @@ -147,7 +147,7 @@ function FullTimeStep(::SimulationMetaData{D,T,S,K,B,L}, SimKernel, SimConstants @inbounds for i in eachindex(Velocityₙ⁺) max_velocity = max(max_velocity, norm(Velocityₙ⁺[i])) end - shift_scale = -SimConstants.CFL * (max_velocity / SimConstants.c₀) * (2 * SimKernel.h)^2 + shift_scale = -SimConstants.ShiftCFL * (max_velocity / SimConstants.c₀) * (2 * SimKernel.h)^2 @inbounds @simd ivdep for i in eachindex(Position) MotionLimiterFactor = MotionLimiterValue(AccelerationScalarType, ParticleType[i]) GravityFactor = GravityFactorValue(AccelerationScalarType, ParticleType[i])