AhmedSalih3d · AhmedSalih3d · Jul 5, 2026 · Jul 5, 2026 · Jul 5, 2026 · Jul 5, 2026
diff --git a/README.md b/README.md
@@ -26,8 +26,9 @@ The project demonstrates how to assemble a small SPH solver with Julia. It focus
 - **Wendland quintic kernel** – simple and stable without tensile corrections.
 - **Symplectic time stepping** – choose between symplectic two-loop and single-loop midpoint updates.
 
-Time-stepping behavior is selected via `RunSimulation(..., SimTimeStepping=...)` with either
-`SymplecticTimeStepping()` or `SingleNeighborTimeStepping()` depending on the desired update path.
+Time-stepping behavior is encoded in the `SimulationMetaData` type. The shorter metadata constructors default to
+`SingleNeighborTimeStepping`; use the seven-parameter form with `SymplecticTimeStepping` when the symplectic
+two-loop update path is desired.
 
 ## Folder Structure
 

diff --git a/example/Dambreak2dMDBC.jl b/example/Dambreak2dMDBC.jl
@@ -79,7 +79,6 @@ let
         SimParticles         = SimParticles,
         SimViscosity         = ArtificialViscosity(),
         SimDensityDiffusion  = LinearDensityDiffusion(),
-        SimTimeStepping      = SingleNeighborTimeStepping(),
         ParticleNormalsPath  = "./input/dam_break_2d/DamBreak2d_Dp0.02_MDBC_GhostNodes_ThreeLayers.csv"
     )
 end
diff --git a/example/Dambreak3d.jl b/example/Dambreak3d.jl
@@ -73,6 +73,5 @@ let
         SimParticles       = SimParticles,
         SimViscosity       = SimViscosity,
         SimDensityDiffusion= SimDensityDiffusion,
-        SimTimeStepping    = SingleNeighborTimeStepping()
     )
 end
diff --git a/example/DucklingMDBC.jl b/example/DucklingMDBC.jl
@@ -60,7 +60,6 @@ let
         SimKernel           = SimKernel,
         SimViscosity        = SimViscosity,
         SimDensityDiffusion = SimDensityDiffusion,
-        SimTimeStepping     = SingleNeighborTimeStepping(),
         ParticleNormalsPath = "./input/case_duckling_mdbc/CaseDuckling_Dp$(SimConstantsWedge.dx)_GhostNodes.csv"
     )
 

diff --git a/example/MovingSquare2d.jl b/example/MovingSquare2d.jl
@@ -82,6 +82,5 @@ let
         SimKernel           = SimKernel,
         SimViscosity        = LaminarSPS(),
         SimDensityDiffusion = LinearDensityDiffusion(),
-        SimTimeStepping     = SingleNeighborTimeStepping()
     )
 end
diff --git a/example/StillWedgeMDBC.jl b/example/StillWedgeMDBC.jl
@@ -58,7 +58,6 @@ let
         SimParticles        = SimParticles,
         SimViscosity        = ArtificialViscosity(),
         SimDensityDiffusion = LinearDensityDiffusion(),
-        SimTimeStepping     = SingleNeighborTimeStepping(),
         ParticleNormalsPath = "./input/still_wedge_mdbc/StillWedge_Dp$(SimConstantsWedge.dx)_GhostNodes_Correct.csv"
     )
 

diff --git a/example/StillWedgeMiddleSquareMDBC.jl b/example/StillWedgeMiddleSquareMDBC.jl
@@ -63,7 +63,6 @@ let
         SimParticles        = SimParticles,
         SimViscosity        = ArtificialViscosity(),
         SimDensityDiffusion = LinearDensityDiffusion(),
-        SimTimeStepping     = SingleNeighborTimeStepping(),
         ParticleNormalsPath = "./input/still_wedge_middle_square_mdbc/StillWedge_MiddleSquare_Dp$(SimConstantsWedge.dx)_GhostNodes.csv"
     )
 

diff --git a/src/SPHCellList.jl b/src/SPHCellList.jl
@@ -688,8 +688,88 @@ using LinearAlgebra
 
     # Per-particle local Δx removed: use single scalar `SimMetaData.Δx`.
 
+    function InitialNeighborLoop!(::SimulationMetaData{D,T,S,K,B,L,SymplecticTimeStepping}, _args...) where {D,T,S<:ShiftingMode,K<:KernelOutputMode,B<:MDBCMode,L<:LogMode}
+        return nothing
+    end
+
+    function InitialNeighborLoop!(SimMetaData::SimulationMetaData{D,T,S,K,B,L,SingleNeighborTimeStepping},
+                                  SimDensityDiffusion, SimViscosity, SimKernel, SimConstants,
+                                  SimParticles, ParticleRanges, CellListIndices, NeighborCellLists,
+                                  dρdtI, ∇Cᵢ, ∇◌rᵢ, AccelerationMax, UniqueCells) where {D,T,S<:ShiftingMode,K<:KernelOutputMode,B<:MDBCMode,L<:LogMode}
+        @timeit SimMetaData.HourGlass "00 Init Pressure"                          Pressure!(SimParticles.Pressure, SimParticles.Density, SimConstants)
+        @timeit SimMetaData.HourGlass "00a Init MDBC"                             ApplyMDBCBeforeHalf!(SimMetaData, SimKernel, SimConstants, SimParticles, ParticleRanges, UniqueCells)
+        @timeit SimMetaData.HourGlass "00b Init NeighborLoop" NeighborLoopPerParticle!(
+            SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData,
+            SimConstants, SimParticles, ParticleRanges, CellListIndices,
+            NeighborCellLists, dρdtI, SimParticles.Acceleration, ∇Cᵢ, ∇◌rᵢ, AccelerationMax,
+        )
+
+        return nothing
+    end
+
+    function HalfStepNeighborLoops!(SimMetaData::SimulationMetaData{D,T,S,K,B,L,SymplecticTimeStepping},
+                                    SimDensityDiffusion, SimViscosity, SimKernel, SimConstants,
+                                    SimParticles, ParticleRanges, UniqueCells, CellListIndices,
+                                    NeighborCellLists, dρdtI, Velocityₙ⁺, Positionₙ⁺, ρₙ⁺,
+                                    ∇Cᵢ, ∇◌rᵢ, AccelerationMax, dt₂, ParticleType,
+                                    MotionDefinition) where {D,T,S<:ShiftingMode,K<:KernelOutputMode,B<:MDBCMode,L<:LogMode}
+        @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 "05 Update To Half TimeStep"               HalfTimeStep(SimMetaData, SimConstants, SimParticles, Positionₙ⁺, Velocityₙ⁺, ρₙ⁺, dρdtI, dt₂)
+
+        @timeit SimMetaData.HourGlass "06 Half LimitDensityAtBoundary"           LimitDensityAtBoundary!(ρₙ⁺, SimConstants.ρ₀, ParticleType)
+
+        @timeit SimMetaData.HourGlass "Motion"                                   ProgressMotion(SimParticles, dt₂, MotionDefinition, SimMetaData)
+
+        @timeit SimMetaData.HourGlass "07 Pressure"                              Pressure!(SimParticles.Pressure, ρₙ⁺, SimConstants)
+        @timeit SimMetaData.HourGlass "08 Second NeighborLoop" NeighborLoopPerParticle!(
+            SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData,
+            SimConstants, SimParticles, ParticleRanges, CellListIndices,
+            NeighborCellLists, dρdtI, SimParticles.Acceleration, ∇Cᵢ, ∇◌rᵢ, AccelerationMax,
+            Position = Positionₙ⁺,
+            Density = ρₙ⁺,
+            Velocity = Velocityₙ⁺,
+        )
+
+        return nothing
+    end
+
+    function HalfStepNeighborLoops!(SimMetaData::SimulationMetaData{D,T,S,K,B,L,SingleNeighborTimeStepping},
+                                    SimDensityDiffusion, SimViscosity, SimKernel, SimConstants,
+                                    SimParticles, ParticleRanges, UniqueCells, CellListIndices,
+                                    NeighborCellLists, dρdtI, Velocityₙ⁺, Positionₙ⁺, ρₙ⁺,
+                                    ∇Cᵢ, ∇◌rᵢ, AccelerationMax, dt₂, ParticleType,
+                                    MotionDefinition) where {D,T,S<:ShiftingMode,K<:KernelOutputMode,B<:MDBCMode,L<:LogMode}
+        @timeit SimMetaData.HourGlass "02 Apply MDBC before Half TimeStep"       ApplyMDBCBeforeHalf!(SimMetaData, SimKernel, SimConstants, SimParticles, ParticleRanges, UniqueCells)
+
+        @timeit SimMetaData.HourGlass "03 Update To Half TimeStep"               HalfTimeStep(SimMetaData, SimConstants, SimParticles, Positionₙ⁺, Velocityₙ⁺, ρₙ⁺, dρdtI, dt₂)
+
+        @timeit SimMetaData.HourGlass "04 Half LimitDensityAtBoundary"           LimitDensityAtBoundary!(ρₙ⁺, SimConstants.ρ₀, ParticleType)
+
+        @timeit SimMetaData.HourGlass "Motion"                                   ProgressMotion(SimParticles, dt₂, MotionDefinition, SimMetaData)
+
+        @timeit SimMetaData.HourGlass "05 Pressure"                              Pressure!(SimParticles.Pressure, ρₙ⁺, SimConstants)
+        @timeit SimMetaData.HourGlass "06 NeighborLoop" NeighborLoopPerParticle!(
+            SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData,
+            SimConstants, SimParticles, ParticleRanges, CellListIndices,
+            NeighborCellLists, dρdtI, SimParticles.Acceleration, ∇Cᵢ, ∇◌rᵢ, AccelerationMax,
+            Position = Positionₙ⁺,
+            Density = ρₙ⁺,
+            Velocity = Velocityₙ⁺,
+        )
+
+        return nothing
+    end
+
     @inbounds function SimulationLoop(SimDensityDiffusion::SDD, SimViscosity::SV, SimKernel,
-                                      SimMetaData::SimulationMetaData{Dimensions, FloatType, SMode, KMode, BMode, LMode},
+                                      SimMetaData::SimulationMetaData{Dimensions, FloatType, SMode, KMode, BMode, LMode, TMode},
                                       SimConstants, SimParticles, FullStencil,
                                       ParticleRanges, UniqueCells, CellListIndices,
                                       SortingScratchSpace,
@@ -706,6 +786,7 @@ using LinearAlgebra
                                       }) where {
                                                 Dimensions, FloatType, SMode, KMode,
                                                 BMode, LMode,
+                                                TMode<:TimeSteppingMode,
                                                 SDD<:SPHDensityDiffusion,
                                                 SV<:SPHViscosity}
         ParticleType   = SimParticles.Type
@@ -717,8 +798,6 @@ using LinearAlgebra
         UniqueCellsView = view(UniqueCells, 1:SimMetaData.IndexCounter)
         # This code here is to initialize the first time step for each simulation loop
         dt = SimConstants.CFL * (SimKernel.h / SimConstants.c₀)
-        TimeSteppingMode = SimMetaData.TimeSteppingMode
-
         @no_escape begin
             AccelerationMax = @alloc(FloatType, length(SimParticles.Position))
             dt₂ = dt * 0.5
@@ -727,15 +806,9 @@ using LinearAlgebra
             UniqueCellsView = view(UniqueCells, 1:SimMetaData.IndexCounter)
             BuildNeighborCellLists!(NeighborCellLists, FullStencil, UniqueCellsView, ParticleRanges)
 
-            if TimeSteppingMode isa SingleNeighborTimeStepping
-                @timeit SimMetaData.HourGlass "00 Init Pressure"                          Pressure!(SimParticles.Pressure, SimParticles.Density, SimConstants)
-                @timeit SimMetaData.HourGlass "00a Init MDBC"                             ApplyMDBCBeforeHalf!(SimMetaData, SimKernel, SimConstants, SimParticles, ParticleRanges, UniqueCells)
-                @timeit SimMetaData.HourGlass "00b Init NeighborLoop" NeighborLoopPerParticle!(
-                    SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData,
-                    SimConstants, SimParticles, ParticleRanges, CellListIndices,
-                    NeighborCellLists, dρdtI, SimParticles.Acceleration, ∇Cᵢ, ∇◌rᵢ, AccelerationMax,
-                )
-            end
+            InitialNeighborLoop!(SimMetaData, SimDensityDiffusion, SimViscosity, SimKernel, SimConstants,
+                                 SimParticles, ParticleRanges, CellListIndices, NeighborCellLists,
+                                 dρdtI, ∇Cᵢ, ∇◌rᵢ, AccelerationMax, UniqueCells)
 
             NextOutputTime = next_output_time(SimMetaData)
             while SimMetaData.TotalTime <= NextOutputTime
@@ -762,50 +835,10 @@ using LinearAlgebra
 
                 @timeit SimMetaData.HourGlass "Motion"                                   ProgressMotion(SimParticles, dt₂, MotionDefinition, SimMetaData)
 
-                if TimeSteppingMode isa SymplecticTimeStepping
-                    @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 "05 Update To Half TimeStep"               HalfTimeStep(SimMetaData, SimConstants, SimParticles, Positionₙ⁺, Velocityₙ⁺, ρₙ⁺, dρdtI, dt₂)
-
-                    @timeit SimMetaData.HourGlass "06 Half LimitDensityAtBoundary"           LimitDensityAtBoundary!(ρₙ⁺, SimConstants.ρ₀, ParticleType)
-
-                    @timeit SimMetaData.HourGlass "Motion"                                   ProgressMotion(SimParticles, dt₂, MotionDefinition, SimMetaData)
-
-                    @timeit SimMetaData.HourGlass "07 Pressure"                              Pressure!(SimParticles.Pressure, ρₙ⁺, SimConstants)
-                    @timeit SimMetaData.HourGlass "08 Second NeighborLoop" NeighborLoopPerParticle!(
-                        SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData,
-                        SimConstants, SimParticles, ParticleRanges, CellListIndices,
-                        NeighborCellLists, dρdtI, SimParticles.Acceleration, ∇Cᵢ, ∇◌rᵢ, AccelerationMax,
-                        Position = Positionₙ⁺,
-                        Density = ρₙ⁺,
-                        Velocity = Velocityₙ⁺,
-                    )
-                else
-                    @timeit SimMetaData.HourGlass "02 Apply MDBC before Half TimeStep"       ApplyMDBCBeforeHalf!(SimMetaData, SimKernel, SimConstants, SimParticles, ParticleRanges, UniqueCells)
-
-                    @timeit SimMetaData.HourGlass "03 Update To Half TimeStep"               HalfTimeStep(SimMetaData, SimConstants, SimParticles, Positionₙ⁺, Velocityₙ⁺, ρₙ⁺, dρdtI, dt₂)
-
-                    @timeit SimMetaData.HourGlass "04 Half LimitDensityAtBoundary"           LimitDensityAtBoundary!(ρₙ⁺, SimConstants.ρ₀, ParticleType)
-
-                    @timeit SimMetaData.HourGlass "Motion"                                   ProgressMotion(SimParticles, dt₂, MotionDefinition, SimMetaData)
-
-                    @timeit SimMetaData.HourGlass "05 Pressure"                              Pressure!(SimParticles.Pressure, ρₙ⁺, SimConstants)
-                    @timeit SimMetaData.HourGlass "06 NeighborLoop" NeighborLoopPerParticle!(
-                        SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData,
-                        SimConstants, SimParticles, ParticleRanges, CellListIndices,
-                        NeighborCellLists, dρdtI, SimParticles.Acceleration, ∇Cᵢ, ∇◌rᵢ, AccelerationMax,
-                        Position = Positionₙ⁺,
-                        Density = ρₙ⁺,
-                        Velocity = Velocityₙ⁺,
-                    )
-                end
+                HalfStepNeighborLoops!(SimMetaData, SimDensityDiffusion, SimViscosity, SimKernel, SimConstants,
+                                       SimParticles, ParticleRanges, UniqueCells, CellListIndices,
+                                       NeighborCellLists, dρdtI, Velocityₙ⁺, Positionₙ⁺, ρₙ⁺,
+                                       ∇Cᵢ, ∇◌rᵢ, AccelerationMax, dt₂, ParticleType, MotionDefinition)
 
                 @timeit SimMetaData.HourGlass "07 Final Density"                         DensityEpsi!(SimParticles.Density, dρdtI, ρₙ⁺, dt)
 
@@ -821,24 +854,20 @@ using LinearAlgebra
 
         return nothing
     end
-
     ###===
     function RunSimulation(;SimGeometry::Vector{Geometry{Dimensions, FloatType}}, #Don't further specify type for now
-        SimMetaData::SimulationMetaData{Dimensions, FloatType, SMode, KMode, BMode, LMode},
+        SimMetaData::SimulationMetaData{Dimensions, FloatType, SMode, KMode, BMode, LMode, TMode},
         SimConstants::SimulationConstants,
         SimKernel::SPHKernelInstance,
         SimLogger::SimulationLogger,
         SimParticles::StructArray,
         SimViscosity::SV,
         SimDensityDiffusion::SDD,
-        SimTimeStepping::TimeSteppingMode,
         ParticleNormalsPath::Union{Nothing,String} = nothing
-        ) where {Dimensions,FloatType,SMode,KMode,BMode,LMode,SV<:SPHViscosity,SDD<:SPHDensityDiffusion}
+        ) where {Dimensions,FloatType,SMode,KMode,BMode,LMode,TMode<:TimeSteppingMode,SV<:SPHViscosity,SDD<:SPHDensityDiffusion}
 
         NumberOfPoints = length(SimParticles)
 
-        SimMetaData.TimeSteppingMode = SimTimeStepping
-
         dρdtI, Velocityₙ⁺, Positionₙ⁺, ρₙ⁺, ∇Cᵢ, ∇◌rᵢ = AllocateSupportDataStructures(SimMetaData, SimParticles.Position)
 
         LoadMDBCNormals!(SimMetaData, SimParticles, ParticleNormalsPath)

diff --git a/src/SimulationMetaDataConfiguration.jl b/src/SimulationMetaDataConfiguration.jl
@@ -35,7 +35,8 @@ struct SingleNeighborTimeStepping <: TimeSteppingMode end
                                            SMode <: ShiftingMode,
                                            KMode <: KernelOutputMode,
                                            BMode <: MDBCMode,
-                                           LMode <: LogMode}
+                                           LMode <: LogMode,
+                                           TMode <: TimeSteppingMode}
     SimulationName::String
     SaveLocation::String
     HourGlass::TimerOutput                  = TimerOutput()
@@ -55,8 +56,10 @@ struct SingleNeighborTimeStepping <: TimeSteppingMode end
     ExportGridCellParticleCounts::Bool      = false
     OpenLogFile::Bool                       = true
     Δx::FloatType                           = zero(FloatType)
-    TimeSteppingMode::TimeSteppingMode      = SingleNeighborTimeStepping()
 end
+
+SimulationMetaData{D,T,S,K,B,L}(; kwargs...) where {D,T,S<:ShiftingMode,K<:KernelOutputMode,B<:MDBCMode,L<:LogMode} =
+    SimulationMetaData{D,T,S,K,B,L,SingleNeighborTimeStepping}(; kwargs...)
 SimulationMetaData{D,T,S,K,B}(; kwargs...) where {D,T,S<:ShiftingMode,K<:KernelOutputMode,B<:MDBCMode} =
     SimulationMetaData{D,T,S,K,B,NoLog}(; kwargs...)
 SimulationMetaData{D,T,S,K}(; kwargs...) where {D,T,S<:ShiftingMode,K<:KernelOutputMode} =