diff --git a/README.md b/README.md index 4d4dc776..7ac7d317 100644 --- a/README.md +++ b/README.md @@ -88,6 +88,9 @@ 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` 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 `ParticleNeighborsPerCell` array that includes each cell's particle count minus diff --git a/src/PreProcess.jl b/src/PreProcess.jl index 1e7630ca..4a999ad7 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 1cb20127..ab5712db 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 @@ -379,6 +385,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 +402,8 @@ using LinearAlgebra vⱼ = Velocity[j] vᵢⱼ = vᵢ - vⱼ density_symmetric_term = dot(-vᵢⱼ, ∇ᵢWᵢⱼ) - dρdt⁺ = -ρᵢ * (m₀ / ρⱼ) * density_symmetric_term + 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) @@ -405,7 +413,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 +434,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 +451,8 @@ using LinearAlgebra vⱼ = Velocity[j] vᵢⱼ = vᵢ - vⱼ density_symmetric_term = dot(-vᵢⱼ, ∇ᵢWᵢⱼ) - dρdt⁺ = -ρᵢ * (m₀ / ρⱼ) * density_symmetric_term + 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) @@ -452,7 +462,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 +485,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 +502,8 @@ using LinearAlgebra vⱼ = Velocity[j] vᵢⱼ = vᵢ - vⱼ density_symmetric_term = dot(-vᵢⱼ, ∇ᵢWᵢⱼ) - dρdt⁺ = -ρᵢ * (m₀ / ρⱼ) * density_symmetric_term + 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) @@ -501,7 +513,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 +523,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 +541,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 +558,8 @@ using LinearAlgebra vⱼ = Velocity[j] vᵢⱼ = vᵢ - vⱼ density_symmetric_term = dot(-vᵢⱼ, ∇ᵢWᵢⱼ) - dρdt⁺ = -ρᵢ * (m₀ / ρⱼ) * density_symmetric_term + 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) @@ -555,7 +569,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 +577,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 @@ -578,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 @@ -610,11 +626,14 @@ using LinearAlgebra first_column..., ((xⱼᵢ * first_column')')... ) + + kernel_sum = VⱼWᵢⱼ + div_pos = Vⱼ * dot(Position[j] - GhostPoints[i], ∇ᵢWᵢⱼ) end end - return bΔ, AΔ + 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} @@ -629,36 +648,109 @@ using LinearAlgebra 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)) + DivPos = @alloc(T, length(SimParticles.Position)) UniqueCellsView = view(UniqueCells, 1:SimMetaData.IndexCounter) - NeighborLoopMDBC!(SimKernel, SimMetaData, SimConstants, ParticleRanges, UniqueCellsView, SimParticles, bᵧ, Aᵧ) - ApplyMDBCCorrection(SimConstants, SimParticles, bᵧ, Aᵧ) + 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) + 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) 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 + 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)) + + kernel_sum = KernelSums[i] + ghost_density = ρ₀ + grad_density = zero(SVector{length(Position[i]), eltype(ρ₀)}) + + use_shepherd = kernel_sum < kernel_sum_threshold + use_matrix = false + if !use_shepherd + det_value = det(A) + if abs(det_value) >= det_threshold + condition_number = cond(A) + use_matrix = condition_number < cond_threshold end + use_shepherd = !use_matrix + end + + if use_matrix + ghost_state = A \ bᵧ[i] + 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 && kernel_sum > zero(eltype(kernel_sum)) + ghost_density = first(bᵧ[i]) / kernel_sum + end + + if isnan(ghost_density) + ghost_density = ρ₀ + end + diff = Position[i] - GhostPoints[i] + boundary_density = ghost_density + dot(diff, grad_density) + if isnan(boundary_density) + boundary_density = ρ₀ + end + + 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 + if !isnan(boundary_pressure) + 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 + if isnan(Pressure[i]) + Pressure[i] = EquationOfStateGamma7(Density[i], c₀, ρ₀) end end end diff --git a/src/SPHDensityDiffusionModels.jl b/src/SPHDensityDiffusionModels.jl index 71651654..d0c1ba8f 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ⱼ = ParticleType[j] == Fluid ? 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ⱼ = ParticleType[j] == Fluid ? 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ⱼ = 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 @@ -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..37d02097 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ⱼ = SimParticles.Type[j] == Fluid ? 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ⱼ = SimParticles.Type[j] == Fluid ? 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ⱼ = SimParticles.Type[j] == Fluid ? 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ⱼ