Skip to content

Add SIMD helpers to SPHCellList #169

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
AhmedSalih3d wants to merge 2 commits into
base: codex/refactor-to-avoid-if-statement
Choose a base branch
from
Open

Add SIMD helpers to SPHCellList #169

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
658c4a3
Add SIMD helpers to cell list
AhmedSalih3d Jan 23, 2026
5409ece
Expand SIMD dot usage
AhmedSalih3d Jan 23, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 2 additions & 0 deletions Project.toml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -16,6 +16,7 @@ Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
LoggingExtras = "e6f89c97-d47a-5376-807f-9c37f3926c36"
Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a"
Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
SIMD = "fdea26ae-647d-5447-a871-4b548cad5224"
StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
StructArrays = "09ab397b-f2b6-538f-b94a-2f83cf4a842a"
TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
Expand All @@ -31,6 +32,7 @@ InteractiveUtils = "1.11.0"
LibGit2 = "1.11.0"
LoggingExtras = "^1.2.0"
Parameters = "^0.12.3"
SIMD = "3.7.2"
StaticArrays = "^1.9.15"
StructArrays = "^0.7.2"
TimerOutputs = "^0.5.29"
Expand Down
26 changes: 25 additions & 1 deletion src/AuxiliaryFunctions.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2,8 +2,10 @@ module AuxiliaryFunctions
using StaticArrays
using Base.Threads
using HDF5
using SIMD: Vec
import LinearAlgebra: dot

export ResetArrays!, to_3d, CloseHDFVTKManually, CleanUpSimulationFolder
export ResetArrays!, to_3d, CloseHDFVTKManually, CleanUpSimulationFolder, DotSimd

"""
ResetArrays!(arrays...)
Expand All @@ -19,6 +21,28 @@ Convert a vector of 2D `SVector`s to 3D by appending a zero z-component.
"""
@inline to_3d(vec_2d) = [SVector(v..., 0.0) for v in vec_2d]

"""
DotSimd(a, b)

SIMD-accelerated dot product for `SVector` inputs where it is practical,
falling back to `dot` for other sizes and types.
"""
@inline function DotSimd(a::SVector{2,T}, b::SVector{2,T}) where {T<:AbstractFloat}
veca = Vec{2,T}(Tuple(a))
vecb = Vec{2,T}(Tuple(b))
return sum(veca * vecb)
end

@inline function DotSimd(a::SVector{3,T}, b::SVector{3,T}) where {T<:AbstractFloat}
veca = Vec{4,T}(a[1], a[2], a[3], zero(T))
vecb = Vec{4,T}(b[1], b[2], b[3], zero(T))
return sum(veca * vecb)
end

@inline function DotSimd(a::SVector{D,T}, b::SVector{D,T}) where {D,T}
return dot(a, b)
end

"""
to_3d!(dest, src)

Expand Down
22 changes: 11 additions & 11 deletions src/SPHCellList.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -381,7 +381,7 @@ using LinearAlgebra
@unpack h⁻¹, H², h = SimKernel

xᵢⱼ = Position[i] - Position[j]
xᵢⱼ² = dot(xᵢⱼ, xᵢⱼ)
xᵢⱼ² = DotSimd(xᵢⱼ, xᵢⱼ)
if xᵢⱼ² <= H²
dᵢⱼ = sqrt(abs(xᵢⱼ²))
dᵢⱼ² = dᵢⱼ^2
Expand All @@ -394,7 +394,7 @@ using LinearAlgebra
vᵢ = Velocity[i]
vⱼ = Velocity[j]
vᵢⱼ = vᵢ - vⱼ
density_symmetric_term = dot(-vᵢⱼ, ∇ᵢWᵢⱼ)
density_symmetric_term = DotSimd(-vᵢⱼ, ∇ᵢWᵢⱼ)
dρdt⁺ = -ρᵢ * (m₀ / ρⱼ) * density_symmetric_term

Dᵢ, _ = compute_density_diffusion(SimDensityDiffusion, SimKernel, SimConstants, SimParticles, xᵢⱼ, ∇ᵢWᵢⱼ, dᵢⱼ², i, j, ParticleType)
Expand Down Expand Up @@ -428,7 +428,7 @@ using LinearAlgebra
@unpack h⁻¹, H², h = SimKernel

xᵢⱼ = Position[i] - Position[j]
xᵢⱼ² = dot(xᵢⱼ, xᵢⱼ)
xᵢⱼ² = DotSimd(xᵢⱼ, xᵢⱼ)
if xᵢⱼ² <= H²
dᵢⱼ = sqrt(abs(xᵢⱼ²))
dᵢⱼ² = dᵢⱼ^2
Expand All @@ -441,7 +441,7 @@ using LinearAlgebra
vᵢ = Velocity[i]
vⱼ = Velocity[j]
vᵢⱼ = vᵢ - vⱼ
density_symmetric_term = dot(-vᵢⱼ, ∇ᵢWᵢⱼ)
density_symmetric_term = DotSimd(-vᵢⱼ, ∇ᵢWᵢⱼ)
dρdt⁺ = -ρᵢ * (m₀ / ρⱼ) * density_symmetric_term

Dᵢ, _ = compute_density_diffusion(SimDensityDiffusion, SimKernel, SimConstants, SimParticles, xᵢⱼ, ∇ᵢWᵢⱼ, dᵢⱼ², i, j, ParticleType)
Expand Down Expand Up @@ -477,7 +477,7 @@ using LinearAlgebra
@unpack h⁻¹, H², h = SimKernel

xᵢⱼ = Position[i] - Position[j]
xᵢⱼ² = dot(xᵢⱼ, xᵢⱼ)
xᵢⱼ² = DotSimd(xᵢⱼ, xᵢⱼ)
if xᵢⱼ² <= H²
dᵢⱼ = sqrt(abs(xᵢⱼ²))
dᵢⱼ² = dᵢⱼ^2
Expand All @@ -490,7 +490,7 @@ using LinearAlgebra
vᵢ = Velocity[i]
vⱼ = Velocity[j]
vᵢⱼ = vᵢ - vⱼ
density_symmetric_term = dot(-vᵢⱼ, ∇ᵢWᵢⱼ)
density_symmetric_term = DotSimd(-vᵢⱼ, ∇ᵢWᵢⱼ)
dρdt⁺ = -ρᵢ * (m₀ / ρⱼ) * density_symmetric_term

Dᵢ, _ = compute_density_diffusion(SimDensityDiffusion, SimKernel, SimConstants, SimParticles, xᵢⱼ, ∇ᵢWᵢⱼ, dᵢⱼ², i, j, ParticleType)
Expand All @@ -511,7 +511,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₀ / ρⱼ) * DotSimd(-xᵢⱼ, ∇ᵢWᵢⱼ) * MotionLimiterCondition
end

return dρdt_acc, acc_acc, kernel_acc, kernel_grad_acc, shift_c_acc, shift_r_acc
Expand All @@ -531,7 +531,7 @@ using LinearAlgebra
@unpack h⁻¹, H², h = SimKernel

xᵢⱼ = Position[i] - Position[j]
xᵢⱼ² = dot(xᵢⱼ, xᵢⱼ)
xᵢⱼ² = DotSimd(xᵢⱼ, xᵢⱼ)
if xᵢⱼ² <= H²
dᵢⱼ = sqrt(abs(xᵢⱼ²))
dᵢⱼ² = dᵢⱼ^2
Expand All @@ -544,7 +544,7 @@ using LinearAlgebra
vᵢ = Velocity[i]
vⱼ = Velocity[j]
vᵢⱼ = vᵢ - vⱼ
density_symmetric_term = dot(-vᵢⱼ, ∇ᵢWᵢⱼ)
density_symmetric_term = DotSimd(-vᵢⱼ, ∇ᵢWᵢⱼ)
dρdt⁺ = -ρᵢ * (m₀ / ρⱼ) * density_symmetric_term

Dᵢ, _ = compute_density_diffusion(SimDensityDiffusion, SimKernel, SimConstants, SimParticles, xᵢⱼ, ∇ᵢWᵢⱼ, dᵢⱼ², i, j, ParticleType)
Expand All @@ -563,7 +563,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₀ / ρⱼ) * DotSimd(-xᵢⱼ, ∇ᵢWᵢⱼ) * MotionLimiterCondition
end

return dρdt_acc, acc_acc, shift_c_acc, shift_r_acc
Expand All @@ -585,7 +585,7 @@ using LinearAlgebra

xᵢⱼ = GhostPoints[i] - Position[j]

xᵢⱼ² = dot(xᵢⱼ, xᵢⱼ)
xᵢⱼ² = DotSimd(xᵢⱼ, xᵢⱼ)
if xᵢⱼ² <= H²
dᵢⱼ = sqrt(abs(xᵢⱼ²))
q = clamp(dᵢⱼ * h⁻¹, 0.0, 2.0)
Expand Down
7 changes: 4 additions & 3 deletions src/SPHDensityDiffusionModels.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3,6 +3,7 @@ module SPHDensityDiffusionModels
using StaticArrays, LinearAlgebra, Parameters
using ..SimulationEquations
using ..SimulationGeometry
using ..AuxiliaryFunctions
#---------------------------------------------------------------
# Exported
#---------------------------------------------------------------
Expand Down Expand Up @@ -80,7 +81,7 @@ struct ZeroGravityLinearDensityDiffusion <: SPHDensityDiffusion end
ρⱼᵢ = ρⱼ - ρᵢ
ψᵢⱼ = 2 * ρⱼᵢ * (-xᵢⱼ) * invdᵢⱼ²η²

Dᵢ = δᵩ * h * c₀ * (m₀/ρⱼ) * dot(ψᵢⱼ, ∇ᵢWᵢⱼ)
Dᵢ = δᵩ * h * c₀ * (m₀/ρⱼ) * DotSimd(ψᵢⱼ, ∇ᵢWᵢⱼ)
Dⱼ = -Dᵢ


Expand Down Expand Up @@ -130,7 +131,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₀/ρⱼ) * DotSimd(ψᵢⱼ, ∇ᵢWᵢⱼ) * MotionLimiterCondition
Dⱼ = -Dᵢ

return Dᵢ, Dⱼ
Expand Down Expand Up @@ -182,7 +183,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₀/ρⱼ) * DotSimd(ψᵢⱼ, ∇ᵢWᵢⱼ) * MotionLimiterCondition
Dⱼ = -Dᵢ

return Dᵢ, Dⱼ
Expand Down
5 changes: 3 additions & 2 deletions src/SPHViscosityModels.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,6 +1,7 @@
module SPHViscosityModels

using StaticArrays, LinearAlgebra, Parameters
using ..AuxiliaryFunctions

export SPHViscosity, ZeroViscosity, ArtificialViscosity, Laminar, LaminarSPS, compute_viscosity

Expand Down Expand Up @@ -61,7 +62,7 @@ end
ρᵢ = SimParticles.Density[i]
ρⱼ = SimParticles.Density[j]

v_dot_x = dot(vᵢⱼ, xᵢⱼ)
v_dot_x = DotSimd(vᵢⱼ, xᵢⱼ)
if v_dot_x < 0
ρ̄ = 0.5 * (ρᵢ + ρⱼ)
μᵢⱼ = h * v_dot_x / (d² + η²)
Expand All @@ -82,7 +83,7 @@ end
ρᵢ = SimParticles.Density[i]
ρⱼ = SimParticles.Density[j]

term = (4 * m₀ * ν₀ * dot(xᵢⱼ, ∇ᵢWᵢⱼ)) / ((ρᵢ + ρⱼ) + (d² + η²))
term = (4 * m₀ * ν₀ * DotSimd(xᵢⱼ, ∇ᵢWᵢⱼ)) / ((ρᵢ + ρⱼ) + (d² + η²))
return term * vᵢⱼ, -term * vᵢⱼ
end

Expand Down