diff --git a/README.md b/README.md index 4d4dc776..09984c5b 100644 --- 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 index 94a5d4b8..2441dc2d 100644 --- 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 index b6ca0ee5..157af341 100644 --- 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 index d7ddf44c..7d6f7444 100644 --- 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 index 71904f32..64a4101a 100644 --- 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 index 2a81d6f1..14a3abc0 100644 --- 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 index e8fa6c1a..69437f63 100644 --- 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 index 1fb92b2a..d38cc986 100644 --- 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 index b92381bf..bf78fe39 100644 --- 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} =