From cd9e1cdc5ec3c132dc2b93f5148394e552c7bd81 Mon Sep 17 00:00:00 2001
From: AhmedSalih3d <36305327+AhmedSalih3d@users.noreply.github.com>
Date: Fri, 23 Jan 2026 11:04:48 +0100
Subject: [PATCH] Add SIMD density update

---
 Project.toml               |  2 ++
 src/SimulationEquations.jl | 29 ++++++++++++++++++++++++++---
 2 files changed, 28 insertions(+), 3 deletions(-)

diff --git a/Project.toml b/Project.toml
index 258d0262..fbcbf516 100644
--- a/Project.toml
+++ b/Project.toml
@@ -16,6 +16,7 @@ LoggingExtras = "e6f89c97-d47a-5376-807f-9c37f3926c36"
 Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca"
+SIMD = "fdea26ae-647d-5447-a871-4b548cad5224"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 StructArrays = "09ab397b-f2b6-538f-b94a-2f83cf4a842a"
 TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
@@ -35,6 +36,7 @@ LoggingExtras = "1.1.0"
 Parameters = "0.12.3"
 Printf = "1.10"
 ProgressMeter = "1.11.0"
+SIMD = "3.5.0"
 StaticArrays = "1.9.15"
 StructArrays = "0.7.1"
 TimerOutputs = "0.5.29"
diff --git a/src/SimulationEquations.jl b/src/SimulationEquations.jl
index fead76e3..bb50e1a8 100644
--- a/src/SimulationEquations.jl
+++ b/src/SimulationEquations.jl
@@ -5,6 +5,7 @@ export EquationOfState, EquationOfStateGamma7, Pressure!, DensityEpsi!, LimitDen
 using StaticArrays
 using Parameters
 using FastPow
+using SIMD
 
 @inline function EquationOfStateGamma7(ρ,c₀,ρ₀)
     return @fastpow ((c₀^2*ρ₀)/7) * ((ρ/ρ₀)^7 - 1)
@@ -26,10 +27,32 @@ end
 # This is to handle the special factor multiplied on density in the time stepping procedure, when
 # using symplectic time stepping
 @inline function DensityEpsi!(Density, dρdtIₙ⁺,ρₙ⁺,Δt)
-    @inbounds for i in eachindex(Density)
-        epsi = - (dρdtIₙ⁺[i] / ρₙ⁺[i]) * Δt
-        Density[i] *= (2 - epsi) / (2 + epsi)
+    _density_epsi_simd!(Density, dρdtIₙ⁺, ρₙ⁺, Δt)
+    return Density
+end
+
+@inline function _density_epsi_simd!(Density, dρdtIₙ⁺, ρₙ⁺, Δt)
+    T = eltype(Density)
+    lanes = T === Float32 ? 8 : 4
+    W = Vec{lanes, T}
+    n = length(Density)
+    i = 1
+
+    @inbounds while i <= n - lanes + 1
+        density_vec = vload(W, Density, i)
+        dρdt_vec = vload(W, dρdtIₙ⁺, i)
+        ρn_vec = vload(W, ρₙ⁺, i)
+        epsi_vec = -(dρdt_vec / ρn_vec) * T(Δt)
+        ratio_vec = (T(2) - epsi_vec) / (T(2) + epsi_vec)
+        vstore!(Density, i, density_vec * ratio_vec)
+        i += lanes
+    end
+
+    @inbounds for j in i:n
+        epsi = - (dρdtIₙ⁺[j] / ρₙ⁺[j]) * Δt
+        Density[j] *= (2 - epsi) / (2 + epsi)
     end
+    return Density
 end
 
 # This version of the function using !Bool(MotionLimiter) instead of BoundaryBool