Add WCSPH shifting corrections (density/velocity gradient accumulation and application)

src/SPHCellList.jl

@@ -510,8 +510,11 @@ 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
+            kernel_value = @fastpow Wᵢⱼ(SimKernel, q)
+            w_ref = @fastpow Wᵢⱼ(SimKernel, dx * h⁻¹)
+            correction = one(T) + T(0.2) * (kernel_value / w_ref)^T(4)
+            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
@@ -562,8 +565,11 @@ 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
+            kernel_value = @fastpow Wᵢⱼ(SimKernel, q)
+            w_ref = @fastpow Wᵢⱼ(SimKernel, dx * h⁻¹)
+            correction = one(T) + T(0.2) * (kernel_value / w_ref)^T(4)
+            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

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"

src/TimeStepping.jl

@@ -143,21 +143,18 @@ function FullTimeStep(::SimulationMetaData{D,T,S,K,B,L}, SimKernel, SimConstants
     @unpack Position, Velocity, Acceleration = SimParticles
     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
+    max_velocity = zero(eltype(eltype(Velocityₙ⁺)))
+    @inbounds for i in eachindex(Velocityₙ⁺)
+        max_velocity = max(max_velocity, norm(Velocityₙ⁺[i]))
+    end
+    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])
         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))
-        else
-            δxᵢ = -A_FSC * A * SimKernel.h * norm(Velocityₙ⁺[i]) * dt * ∇Cᵢ[i]
-        end
+        δxᵢ = shift_scale * ∇Cᵢ[i]
 
         Position[i]           += (Velocityₙ⁺[i] * dt + δxᵢ) * MotionLimiterFactor
     end