diff --git a/Project.toml b/Project.toml index 3eae26b8..fd2b9e39 100644 --- a/Project.toml +++ b/Project.toml @@ -5,6 +5,7 @@ version = "0.7.0" [deps] Bumper = "8ce10254-0962-460f-a3d8-1f77fea1446e" CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b" +CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" Dates = "ade2ca70-3891-5945-98fb-dc099432e06a" FastPow = "c0e83750-1142-43a8-81cf-6c956b72b4d1" Format = "1fa38f19-a742-5d3f-a2b9-30dd87b9d5f8" @@ -24,6 +25,7 @@ UnicodePlots = "b8865327-cd53-5732-bb35-84acbb429228" [compat] Bumper = "^0.7.1" CSV = "^0.10.15" +CUDA = "^5.4.0" FastPow = "^0.1.0" Format = "^1.3.7" HDF5 = "^0.17.2" diff --git a/src/SPHCUDA.jl b/src/SPHCUDA.jl index 6b9b31bb..9def9b85 100644 --- a/src/SPHCUDA.jl +++ b/src/SPHCUDA.jl @@ -6,6 +6,7 @@ using CUDA using LinearAlgebra using Parameters using StaticArrays +using TimerOutputs using ..SPHKernels using ..SPHViscosityModels @@ -26,6 +27,152 @@ using StructArrays: StructArray CUDAAvailable() = CUDA.functional() +struct LaunchConfig + Threads::Int + Blocks::Int +end + +@inline function KernelLaunchConfig(length) + if length == 0 + return LaunchConfig(1, 1) + end + threads = min(256, length) + return LaunchConfig(threads, cld(length, threads)) +end + +struct CUDAParticleBuffers{D, T} + Position::CuArray{SVector{D, T}, 1} + Density::CuArray{T, 1} + Pressure::CuArray{T, 1} + Velocity::CuArray{SVector{D, T}, 1} + Acceleration::CuArray{SVector{D, T}, 1} + MotionLimiter::CuArray{T, 1} + GravityFactor::CuArray{T, 1} +end + +struct CUDASupportBuffers{D, T} + DρdtI::CuArray{T, 1} + Velocityₙ⁺::CuArray{SVector{D, T}, 1} + Positionₙ⁺::CuArray{SVector{D, T}, 1} + ρₙ⁺::CuArray{T, 1} + ∇Cᵢ::CuArray{SVector{D, T}, 1} + ∇◌rᵢ::CuArray{T, 1} + ΔtViscous::CuArray{T, 1} + ΔtForce::CuArray{T, 1} + ΔxScratch::CuArray{T, 1} +end + +mutable struct CUDAGridBuffers{D} + CellCoords::CuArray{SVector{D, Int}, 1} + BucketCounts::CuArray{Int, 1} + BucketEnds::CuArray{Int, 1} + BucketStarts::CuArray{Int, 1} + BucketWrite::CuArray{Int, 1} + BucketParticles::CuArray{Int, 1} + NeighborOffsets::CuArray{SVector{D, Int}, 1} + BucketCount::Int +end + +struct CUDAMotionBuffers{D, T} + Velocity::CuArray{T, 1} + StartTime::CuArray{T, 1} + Duration::CuArray{T, 1} + Direction::CuArray{SVector{D, T}, 1} + Active::CuArray{Bool, 1} +end + +function BuildCUDAParticleBuffers(SimParticles) + return CUDAParticleBuffers( + CuArray(SimParticles.Position), + CuArray(SimParticles.Density), + CuArray(SimParticles.Pressure), + CuArray(SimParticles.Velocity), + CuArray(SimParticles.Acceleration), + CuArray(SimParticles.MotionLimiter), + CuArray(SimParticles.GravityFactor), + ) +end + +function BuildCUDASupportBuffers(dρdtI, Velocityₙ⁺, Positionₙ⁺, ρₙ⁺, ∇Cᵢ, ∇◌rᵢ) + return CUDASupportBuffers( + CuArray(dρdtI), + CuArray(Velocityₙ⁺), + CuArray(Positionₙ⁺), + CuArray(ρₙ⁺), + CuArray(∇Cᵢ), + CuArray(∇◌rᵢ), + CuArray(similar(dρdtI)), + CuArray(similar(dρdtI)), + CuArray(similar(dρdtI)), + ) +end + +function BuildCUDAMotionBuffers(SimParticles, MotionDefinition, ::Val{Dimensions}, ::Type{FloatType}) where {Dimensions, FloatType} + NumberOfPoints = length(SimParticles.Position) + ZeroVector = zero(eltype(SimParticles.Position)) + MotionVelocity = zeros(FloatType, NumberOfPoints) + MotionStartTime = zeros(FloatType, NumberOfPoints) + MotionDuration = zeros(FloatType, NumberOfPoints) + MotionDirection = Vector{SVector{Dimensions, FloatType}}(undef, NumberOfPoints) + MotionActive = falses(NumberOfPoints) + + if MotionDefinition !== nothing + @inbounds for i in 1:NumberOfPoints + MotionDirection[i] = ZeroVector + if SimParticles.Type[i] == Moving + motion = MotionDefinition[SimParticles.GroupMarker[i]] + if motion !== nothing + MotionVelocity[i] = motion.Velocity + MotionStartTime[i] = motion.StartTime + MotionDuration[i] = motion.Duration + MotionDirection[i] = motion.Direction + MotionActive[i] = true + end + end + end + else + fill!(MotionDirection, ZeroVector) + end + + return CUDAMotionBuffers( + CuArray(MotionVelocity), + CuArray(MotionStartTime), + CuArray(MotionDuration), + CuArray(MotionDirection), + CuArray(MotionActive), + ) +end + +function BuildCUDAGridBuffers(SimParticles, FullStencil, ::Val{D}) where {D} + NumberOfPoints = length(SimParticles.Position) + BucketCount = max(1, 2 * NumberOfPoints) + NeighborOffsets = [SVector{D, Int}(Tuple(offset)) for offset in FullStencil] + return CUDAGridBuffers( + CuArray(similar(SimParticles.Position, SVector{D, Int})), + CuArray(zeros(Int, BucketCount)), + CuArray(zeros(Int, BucketCount)), + CuArray(zeros(Int, BucketCount)), + CuArray(zeros(Int, BucketCount)), + CuArray(zeros(Int, NumberOfPoints)), + CuArray(NeighborOffsets), + BucketCount, + ) +end + +function SyncParticlesToHost!(SimParticles, CUDABuffers::CUDAParticleBuffers) + CUDA.copyto!(SimParticles.Position, CUDABuffers.Position) + CUDA.copyto!(SimParticles.Density, CUDABuffers.Density) + CUDA.copyto!(SimParticles.Pressure, CUDABuffers.Pressure) + CUDA.copyto!(SimParticles.Velocity, CUDABuffers.Velocity) + CUDA.copyto!(SimParticles.Acceleration, CUDABuffers.Acceleration) + return nothing +end + +function SyncPositionsToHost!(SimParticles, CUDABuffers::CUDAParticleBuffers) + CUDA.copyto!(SimParticles.Position, CUDABuffers.Position) + return nothing +end + @inline function BuildNeighborCellRanges(NeighborCellLists, cell_count) starts = zeros(Int, cell_count + 1) total = 0 @@ -46,12 +193,276 @@ CUDAAvailable() = CUDA.functional() return starts, entries end -@inline function BuildCellListIndices(Cells, CellLookup) - indices = similar(Cells, Int) - @inbounds for i in eachindex(Cells) - indices[i] = CellLookupIndex(CellLookup, Cells[i], 1) +function PressureCUDAKernel!(Press, Density, SimConstants) + i = (blockIdx().x - 1) * blockDim().x + threadIdx().x + if i <= length(Press) + ρ = Density[i] + Press[i] = ((SimConstants.c₀^2 * SimConstants.ρ₀) / 7) * ((ρ / SimConstants.ρ₀)^7 - 1) + end + return nothing +end + +function PressureCUDA!(Press, Density, SimConstants, Launch) + @cuda threads=Launch.Threads blocks=Launch.Blocks PressureCUDAKernel!(Press, Density, SimConstants) + return nothing +end + +function DensityEpsiCUDAKernel!(Density, dρdtIₙ⁺, ρₙ⁺, Δt) + i = (blockIdx().x - 1) * blockDim().x + threadIdx().x + if i <= length(Density) + epsi = -(dρdtIₙ⁺[i] / ρₙ⁺[i]) * Δt + Density[i] = Density[i] * (2 - epsi) / (2 + epsi) + end + return nothing +end + +function DensityEpsiCUDA!(Density, dρdtIₙ⁺, ρₙ⁺, Δt, Launch) + @cuda threads=Launch.Threads blocks=Launch.Blocks DensityEpsiCUDAKernel!(Density, dρdtIₙ⁺, ρₙ⁺, Δt) + return nothing +end + +function LimitDensityAtBoundaryCUDAKernel!(Density, ρ₀, MotionLimiter) + i = (blockIdx().x - 1) * blockDim().x + threadIdx().x + if i <= length(Density) + if (Density[i] < ρ₀) * !Bool(MotionLimiter[i]) + Density[i] = ρ₀ + end + end + return nothing +end + +function LimitDensityAtBoundaryCUDA!(Density, ρ₀, MotionLimiter, Launch) + @cuda threads=Launch.Threads blocks=Launch.Blocks LimitDensityAtBoundaryCUDAKernel!(Density, ρ₀, MotionLimiter) + return nothing +end + +function HalfTimeStepCUDAKernel!(Position, Density, Velocity, Acceleration, GravityFactor, MotionLimiter, + Positionₙ⁺, Velocityₙ⁺, ρₙ⁺, dρdtI, dt₂, SimConstants) + i = (blockIdx().x - 1) * blockDim().x + threadIdx().x + if i <= length(Position) + acc = Acceleration[i] + ConstructGravitySVector(Acceleration[i], SimConstants.g * GravityFactor[i]) + Acceleration[i] = acc + limiter = MotionLimiter[i] + Positionₙ⁺[i] = Position[i] + Velocity[i] * dt₂ * limiter + Velocityₙ⁺[i] = Velocity[i] + acc * dt₂ * limiter + ρₙ⁺[i] = Density[i] + dρdtI[i] * dt₂ + end + return nothing +end + +function HalfTimeStepCUDA!(CUDAParticles::CUDAParticleBuffers, CUDASupport::CUDASupportBuffers, dt₂, SimConstants, Launch) + @cuda threads=Launch.Threads blocks=Launch.Blocks HalfTimeStepCUDAKernel!( + CUDAParticles.Position, + CUDAParticles.Density, + CUDAParticles.Velocity, + CUDAParticles.Acceleration, + CUDAParticles.GravityFactor, + CUDAParticles.MotionLimiter, + CUDASupport.Positionₙ⁺, + CUDASupport.Velocityₙ⁺, + CUDASupport.ρₙ⁺, + CUDASupport.DρdtI, + dt₂, + SimConstants, + ) + return nothing +end + +function FullTimeStepCUDAKernel!(Position, Velocity, Acceleration, GravityFactor, MotionLimiter, dt, SimConstants) + i = (blockIdx().x - 1) * blockDim().x + threadIdx().x + if i <= length(Position) + acc = Acceleration[i] + ConstructGravitySVector(Acceleration[i], SimConstants.g * GravityFactor[i]) + Acceleration[i] = acc + limiter = MotionLimiter[i] + Velocity[i] = Velocity[i] + acc * dt * limiter + Position[i] = Position[i] + (((Velocity[i] + (Velocity[i] - acc * dt * limiter)) / 2) * dt) * limiter end - return indices + return nothing +end + +function FullTimeStepCUDA!(CUDAParticles::CUDAParticleBuffers, dt, SimConstants, Launch) + @cuda threads=Launch.Threads blocks=Launch.Blocks FullTimeStepCUDAKernel!( + CUDAParticles.Position, + CUDAParticles.Velocity, + CUDAParticles.Acceleration, + CUDAParticles.GravityFactor, + CUDAParticles.MotionLimiter, + dt, + SimConstants, + ) + return nothing +end + +function ProgressMotionCUDAKernel!(Position, Velocity, MotionVelocity, MotionStartTime, MotionDuration, + MotionDirection, MotionActive, TotalTime, dt₂) + i = (blockIdx().x - 1) * blockDim().x + threadIdx().x + if i <= length(Position) + if MotionActive[i] + should_move = (MotionStartTime[i] <= TotalTime) && (TotalTime <= (MotionStartTime[i] + MotionDuration[i])) + if should_move + Velocity[i] = MotionVelocity[i] * MotionDirection[i] + else + Velocity[i] = zero(Position[i]) + end + Position[i] = Position[i] + Velocity[i] * dt₂ + end + end + return nothing +end + +function ProgressMotionCUDA!(CUDAParticles::CUDAParticleBuffers, MotionBuffers::CUDAMotionBuffers, TotalTime, dt₂, Launch) + @cuda threads=Launch.Threads blocks=Launch.Blocks ProgressMotionCUDAKernel!( + CUDAParticles.Position, + CUDAParticles.Velocity, + MotionBuffers.Velocity, + MotionBuffers.StartTime, + MotionBuffers.Duration, + MotionBuffers.Direction, + MotionBuffers.Active, + TotalTime, + dt₂, + ) + return nothing +end + +function ΔtCUDAKernel!(ViscousBuffer, ForceBuffer, Position, Velocity, Acceleration, SimConstants, SimKernel) + i = (blockIdx().x - 1) * blockDim().x + threadIdx().x + if i <= length(Position) + r = Position[i] + v = Velocity[i] + a = Acceleration[i] + + r_sq = dot(r, r) + ViscousBuffer[i] = abs(SimKernel.h * dot(v, r) / (r_sq + SimKernel.η²)) + + a_sq = dot(a, a) + if a_sq > 0 + a_mag = sqrt(a_sq) + ForceBuffer[i] = sqrt(SimKernel.h / a_mag) + else + ForceBuffer[i] = Inf + end + end + return nothing +end + +function ΔtCUDA(CUDASupport::CUDASupportBuffers, CUDAParticles::CUDAParticleBuffers, SimConstants, SimKernel, Launch) + @cuda threads=Launch.Threads blocks=Launch.Blocks ΔtCUDAKernel!( + CUDASupport.ΔtViscous, + CUDASupport.ΔtForce, + CUDAParticles.Position, + CUDAParticles.Velocity, + CUDAParticles.Acceleration, + SimConstants, + SimKernel, + ) + max_visc = CUDA.reduce(max, CUDASupport.ΔtViscous) + min_force = CUDA.reduce(min, CUDASupport.ΔtForce) + return SimConstants.CFL * min(min_force, SimKernel.h / (SimConstants.c₀ + max_visc)) +end + +function MaxDisplacementCUDAKernel!(Scratch, Positionₙ⁺, Position) + i = (blockIdx().x - 1) * blockDim().x + threadIdx().x + if i <= length(Position) + diff = Positionₙ⁺[i] - Position[i] + Scratch[i] = sqrt(dot(diff, diff)) + end + return nothing +end + +function UpdateΔxCUDA!(Δx, CUDASupport::CUDASupportBuffers, CUDAParticles::CUDAParticleBuffers, Launch) + @cuda threads=Launch.Threads blocks=Launch.Blocks MaxDisplacementCUDAKernel!( + CUDASupport.ΔxScratch, + CUDASupport.Positionₙ⁺, + CUDAParticles.Position, + ) + maxd = CUDA.reduce(max, CUDASupport.ΔxScratch) + return Δx + 4 * maxd +end + +@inline function MapFloorGPU(X, InverseCutOff) + return Int(sign(X)) * unsafe_trunc(Int, muladd(abs(X), InverseCutOff, 0.5)) +end + +@inline function HashCell(Cell) + if length(Cell) == 2 + return Cell[1] * 73856093 ⊻ Cell[2] * 19349663 + end + return Cell[1] * 73856093 ⊻ Cell[2] * 19349663 ⊻ Cell[3] * 83492791 +end + +function ResetBucketsKernel!(BucketCounts, BucketEnds, BucketStarts, BucketWrite) + i = (blockIdx().x - 1) * blockDim().x + threadIdx().x + if i <= length(BucketCounts) + BucketCounts[i] = 0 + BucketEnds[i] = 0 + BucketStarts[i] = 0 + BucketWrite[i] = 0 + end + return nothing +end + +function BuildCellCoordsKernel!(CellCoords, BucketCounts, Position, InverseCutOff, BucketCount) + i = (blockIdx().x - 1) * blockDim().x + threadIdx().x + if i <= length(Position) + pos = Position[i] + CellCoords[i] = SVector{length(pos), Int}(ntuple(d -> MapFloorGPU(pos[d], InverseCutOff), length(pos))) + hash_val = HashCell(CellCoords[i]) + bucket = mod(hash_val, BucketCount) + 1 + CUDA.atomic_add!(BucketCounts, bucket, 1) + end + return nothing +end + +function InitializeBucketStartsKernel!(BucketCounts, BucketEnds, BucketStarts, BucketWrite) + i = (blockIdx().x - 1) * blockDim().x + threadIdx().x + if i <= length(BucketCounts) + count = BucketCounts[i] + BucketStarts[i] = BucketEnds[i] - count + 1 + BucketWrite[i] = BucketStarts[i] + end + return nothing +end + +function FillBucketsKernel!(BucketWrite, BucketParticles, CellCoords, BucketCount) + i = (blockIdx().x - 1) * blockDim().x + threadIdx().x + if i <= length(CellCoords) + hash_val = HashCell(CellCoords[i]) + bucket = mod(hash_val, BucketCount) + 1 + index = CUDA.atomic_add!(BucketWrite, bucket, 1) + BucketParticles[index] = i + end + return nothing +end + +function RebuildNeighborGridCUDA!(CUDAGrid::CUDAGridBuffers, CUDAParticles::CUDAParticleBuffers, SimKernel, Launch) + @cuda threads=Launch.Threads blocks=Launch.Blocks ResetBucketsKernel!( + CUDAGrid.BucketCounts, + CUDAGrid.BucketEnds, + CUDAGrid.BucketStarts, + CUDAGrid.BucketWrite, + ) + @cuda threads=Launch.Threads blocks=Launch.Blocks BuildCellCoordsKernel!( + CUDAGrid.CellCoords, + CUDAGrid.BucketCounts, + CUDAParticles.Position, + SimKernel.H⁻¹, + CUDAGrid.BucketCount, + ) + CUDAGrid.BucketEnds .= CUDA.cumsum(CUDAGrid.BucketCounts) + @cuda threads=Launch.Threads blocks=Launch.Blocks InitializeBucketStartsKernel!( + CUDAGrid.BucketCounts, + CUDAGrid.BucketEnds, + CUDAGrid.BucketStarts, + CUDAGrid.BucketWrite, + ) + @cuda threads=Launch.Threads blocks=Launch.Blocks FillBucketsKernel!( + CUDAGrid.BucketWrite, + CUDAGrid.BucketParticles, + CUDAGrid.CellCoords, + CUDAGrid.BucketCount, + ) + return nothing end @inline function ComputeDensityDiffusionGPU(::ZeroDensityDiffusion, _SimKernel, _SimConstants, @@ -106,92 +517,56 @@ end end function NeighborLoopCUDAKernel!(dρdtI, Acceleration, Position, Density, Pressure, Velocity, - MotionLimiter, CellListIndices, ParticleRanges, ParticleOrder, - NeighborCellStarts, NeighborCellEntries, SimKernel, SimConstants, + MotionLimiter, CellCoords, BucketStarts, BucketEnds, + BucketParticles, NeighborOffsets, BucketCount, SimKernel, SimConstants, SimDensityDiffusion, SimViscosity) i = (blockIdx().x - 1) * blockDim().x + threadIdx().x if i <= length(Position) dρdt_acc = zero(eltype(dρdtI)) acc_acc = zero(Position[i]) - cell_list_index = CellListIndices[i] - same_cell_start = ParticleRanges[cell_list_index] - same_cell_end = ParticleRanges[cell_list_index + 1] - 1 - - @inbounds for j in same_cell_start:same_cell_end - j_index = ParticleOrder[j] - if j_index != i - xᵢⱼ = Position[i] - Position[j_index] - xᵢⱼ² = dot(xᵢⱼ, xᵢⱼ) - if xᵢⱼ² <= SimKernel.H² - dᵢⱼ = sqrt(abs(xᵢⱼ²)) - q = clamp(dᵢⱼ * SimKernel.h⁻¹, zero(eltype(dᵢⱼ)), one(eltype(dᵢⱼ)) * 2) - ∇ᵢWᵢⱼ = ∇Wᵢⱼ(SimKernel, q, xᵢⱼ) - - ρᵢ = Density[i] - ρⱼ = Density[j_index] - vᵢⱼ = Velocity[i] - Velocity[j_index] - - density_symmetric_term = dot(-vᵢⱼ, ∇ᵢWᵢⱼ) - dρdt⁺ = -ρᵢ * (SimConstants.m₀ / ρⱼ) * density_symmetric_term - - Dᵢ = ComputeDensityDiffusionGPU(SimDensityDiffusion, SimKernel, SimConstants, - Density, MotionLimiter, xᵢⱼ, ∇ᵢWᵢⱼ, - dᵢⱼ^2, i, j_index) - - dρdt_acc += dρdt⁺ + Dᵢ - - Pᵢ = Pressure[i] - Pⱼ = Pressure[j_index] - Pfac = (Pᵢ + Pⱼ) / (ρᵢ * ρⱼ) - f_ab = tensile_correction(SimKernel, Pᵢ, ρᵢ, Pⱼ, ρⱼ, q, SimConstants.dx) - dvdt⁺ = -SimConstants.m₀ * (Pfac + f_ab) * ∇ᵢWᵢⱼ - - visc_term = ComputeViscosityGPU(SimViscosity, SimKernel, SimConstants, - Density, xᵢⱼ, vᵢⱼ, ∇ᵢWᵢⱼ, dᵢⱼ^2, i, j_index) - - acc_acc += dvdt⁺ + visc_term - end - end - end - - neighbor_start = NeighborCellStarts[cell_list_index] - neighbor_end = NeighborCellStarts[cell_list_index + 1] - 1 - @inbounds for neighbor_cursor in neighbor_start:neighbor_end - neighbor_index = NeighborCellEntries[neighbor_cursor] - start_index = ParticleRanges[neighbor_index] - end_index = ParticleRanges[neighbor_index + 1] - 1 - for j in start_index:end_index - j_index = ParticleOrder[j] - xᵢⱼ = Position[i] - Position[j_index] - xᵢⱼ² = dot(xᵢⱼ, xᵢⱼ) - if xᵢⱼ² <= SimKernel.H² - dᵢⱼ = sqrt(abs(xᵢⱼ²)) - q = clamp(dᵢⱼ * SimKernel.h⁻¹, zero(eltype(dᵢⱼ)), one(eltype(dᵢⱼ)) * 2) - ∇ᵢWᵢⱼ = ∇Wᵢⱼ(SimKernel, q, xᵢⱼ) - - ρᵢ = Density[i] - ρⱼ = Density[j_index] - vᵢⱼ = Velocity[i] - Velocity[j_index] - - density_symmetric_term = dot(-vᵢⱼ, ∇ᵢWᵢⱼ) - dρdt⁺ = -ρᵢ * (SimConstants.m₀ / ρⱼ) * density_symmetric_term - - Dᵢ = ComputeDensityDiffusionGPU(SimDensityDiffusion, SimKernel, SimConstants, - Density, MotionLimiter, xᵢⱼ, ∇ᵢWᵢⱼ, - dᵢⱼ^2, i, j_index) - - dρdt_acc += dρdt⁺ + Dᵢ - - Pᵢ = Pressure[i] - Pⱼ = Pressure[j_index] - Pfac = (Pᵢ + Pⱼ) / (ρᵢ * ρⱼ) - f_ab = tensile_correction(SimKernel, Pᵢ, ρᵢ, Pⱼ, ρⱼ, q, SimConstants.dx) - dvdt⁺ = -SimConstants.m₀ * (Pfac + f_ab) * ∇ᵢWᵢⱼ - - visc_term = ComputeViscosityGPU(SimViscosity, SimKernel, SimConstants, - Density, xᵢⱼ, vᵢⱼ, ∇ᵢWᵢⱼ, dᵢⱼ^2, i, j_index) - - acc_acc += dvdt⁺ + visc_term + cell = CellCoords[i] + @inbounds for offset_index in eachindex(NeighborOffsets) + neighbor_cell = cell + NeighborOffsets[offset_index] + hash_val = HashCell(neighbor_cell) + bucket = mod(hash_val, BucketCount) + 1 + start_index = BucketStarts[bucket] + end_index = BucketEnds[bucket] + if start_index <= end_index + @inbounds for j in start_index:end_index + j_index = BucketParticles[j] + if CellCoords[j_index] == neighbor_cell && j_index != i + xᵢⱼ = Position[i] - Position[j_index] + xᵢⱼ² = dot(xᵢⱼ, xᵢⱼ) + if xᵢⱼ² <= SimKernel.H² + dᵢⱼ = sqrt(abs(xᵢⱼ²)) + q = clamp(dᵢⱼ * SimKernel.h⁻¹, zero(eltype(dᵢⱼ)), one(eltype(dᵢⱼ)) * 2) + ∇ᵢWᵢⱼ = ∇Wᵢⱼ(SimKernel, q, xᵢⱼ) + + ρᵢ = Density[i] + ρⱼ = Density[j_index] + vᵢⱼ = Velocity[i] - Velocity[j_index] + + density_symmetric_term = dot(-vᵢⱼ, ∇ᵢWᵢⱼ) + dρdt⁺ = -ρᵢ * (SimConstants.m₀ / ρⱼ) * density_symmetric_term + + Dᵢ = ComputeDensityDiffusionGPU(SimDensityDiffusion, SimKernel, SimConstants, + Density, MotionLimiter, xᵢⱼ, ∇ᵢWᵢⱼ, + dᵢⱼ^2, i, j_index) + + dρdt_acc += dρdt⁺ + Dᵢ + + Pᵢ = Pressure[i] + Pⱼ = Pressure[j_index] + Pfac = (Pᵢ + Pⱼ) / (ρᵢ * ρⱼ) + f_ab = tensile_correction(SimKernel, Pᵢ, ρᵢ, Pⱼ, ρⱼ, q, SimConstants.dx) + dvdt⁺ = -SimConstants.m₀ * (Pfac + f_ab) * ∇ᵢWᵢⱼ + + visc_term = ComputeViscosityGPU(SimViscosity, SimKernel, SimConstants, + Density, xᵢⱼ, vᵢⱼ, ∇ᵢWᵢⱼ, dᵢⱼ^2, i, j_index) + + acc_acc += dvdt⁺ + visc_term + end + end end end end @@ -203,20 +578,14 @@ function NeighborLoopCUDAKernel!(dρdtI, Acceleration, Position, Density, Pressu end function NeighborLoopPerParticleCUDA!(SimDensityDiffusion::SDD, SimViscosity::SV, SimKernel, - SimMetaData::SimulationMetaData{D,T,NoShifting,NoKernelOutput,B,L}, - SimConstants, SimParticles, ParticleRanges, - CellLookup, NeighborCellLists, dρdtI, - Acceleration, ∇Cᵢ, - ∇◌rᵢ, max_visc = nothing, - min_dt_force = nothing; - Position = SimParticles.Position, - Density = SimParticles.Density, - Pressure = SimParticles.Pressure, - Velocity = SimParticles.Velocity, - ParticleOrder) where {D,T, - B<:MDBCMode,L<:LogMode, - SDD<:SPHDensityDiffusion, - SV<:SPHViscosity} + SimConstants, CUDAParticles::CUDAParticleBuffers, + CUDAGrid::CUDAGridBuffers, dρdtI, Acceleration, Launch; + Position = CUDAParticles.Position, + Density = CUDAParticles.Density, + Pressure = CUDAParticles.Pressure, + Velocity = CUDAParticles.Velocity) where { + SDD<:SPHDensityDiffusion, + SV<:SPHViscosity} if !CUDAAvailable() error("CUDA is not available; cannot run GPU neighbor loop.") end @@ -228,41 +597,14 @@ function NeighborLoopPerParticleCUDA!(SimDensityDiffusion::SDD, SimViscosity::SV error("CUDA neighbor loop supports ArtificialViscosity or ZeroViscosity.") end - cell_list_indices = BuildCellListIndices(SimParticles.Cells, CellLookup) - cell_count = length(NeighborCellLists) - neighbor_cell_starts, neighbor_cell_entries = BuildNeighborCellRanges(NeighborCellLists, cell_count) - - dρdtI_gpu = CuArray(dρdtI) - acceleration_gpu = CuArray(Acceleration) - position_gpu = CuArray(Position) - density_gpu = CuArray(Density) - pressure_gpu = CuArray(Pressure) - velocity_gpu = CuArray(Velocity) - motion_limiter_gpu = CuArray(SimParticles.MotionLimiter) - cell_list_indices_gpu = CuArray(cell_list_indices) - particle_ranges_gpu = CuArray(ParticleRanges) - particle_order_gpu = CuArray(ParticleOrder) - neighbor_cell_starts_gpu = CuArray(neighbor_cell_starts) - neighbor_cell_entries_gpu = CuArray(neighbor_cell_entries) - - threads = 256 - blocks = cld(length(Position), threads) - @cuda threads=threads blocks=blocks NeighborLoopCUDAKernel!( - dρdtI_gpu, acceleration_gpu, position_gpu, density_gpu, pressure_gpu, - velocity_gpu, motion_limiter_gpu, cell_list_indices_gpu, particle_ranges_gpu, - particle_order_gpu, neighbor_cell_starts_gpu, neighbor_cell_entries_gpu, + @cuda threads=Launch.Threads blocks=Launch.Blocks NeighborLoopCUDAKernel!( + dρdtI, Acceleration, Position, Density, Pressure, + Velocity, CUDAParticles.MotionLimiter, CUDAGrid.CellCoords, + CUDAGrid.BucketStarts, CUDAGrid.BucketEnds, CUDAGrid.BucketParticles, + CUDAGrid.NeighborOffsets, CUDAGrid.BucketCount, SimKernel, SimConstants, SimDensityDiffusion, SimViscosity, ) - CUDA.copyto!(dρdtI, dρdtI_gpu) - CUDA.copyto!(Acceleration, acceleration_gpu) - - if max_visc !== nothing || min_dt_force !== nothing - @inbounds for i in eachindex(Position) - UpdateTimeStepBuffers!(max_visc, min_dt_force, i, Position[i], Velocity[i], Acceleration[i], SimKernel) - end - end - return nothing end @@ -271,99 +613,67 @@ end SimConstants, SimParticles, FullStencil, ParticleRanges, UniqueCells, CellLookup, ParticleOrder, CellOffsets, - NeighborCellLists, dρdtI, Velocityₙ⁺, - Positionₙ⁺, ρₙ⁺, ∇Cᵢ, ∇◌rᵢ, - MotionDefinition::Union{ - Nothing, - AbstractVector{ - Union{ - Nothing, - MotionDetails{Dimensions, FloatType}, - }, - }, - }) where { - Dimensions, FloatType, SMode, KMode, - BMode, LMode, - SDD<:SPHDensityDiffusion, - SV<:SPHViscosity} - @unpack Position, Density, Pressure, Velocity, Acceleration, MotionLimiter, - GroupMarker = SimParticles - ParticleType = SimParticles.Type - ParticleMarker = GroupMarker - GhostPoints = hasproperty(SimParticles, :GhostPoints) ? SimParticles.GhostPoints : nothing - GhostNormals = hasproperty(SimParticles, :GhostNormals) ? SimParticles.GhostNormals : nothing - - UniqueCellsView = view(UniqueCells, 1:SimMetaData.IndexCounter) - dt = Δt(Position, Velocity, Acceleration, SimConstants, SimKernel) - - max_visc = Vector{FloatType}(undef, length(Position)) - min_dt_force = Vector{FloatType}(undef, length(Position)) - + NeighborCellLists, CUDAParticles::CUDAParticleBuffers, + CUDASupport::CUDASupportBuffers, MotionBuffers, + CUDAGrid::CUDAGridBuffers) where { + Dimensions, FloatType, SMode, KMode, + BMode, LMode, + SDD<:SPHDensityDiffusion, + SV<:SPHViscosity} + Launch = KernelLaunchConfig(length(CUDAParticles.Position)) + RebuildNeighborGridCUDA!(CUDAGrid, CUDAParticles, SimKernel, Launch) + dt = ΔtCUDA(CUDASupport, CUDAParticles, SimConstants, SimKernel, Launch) dt₂ = dt * 0.5 while SimMetaData.TotalTime <= next_output_time(SimMetaData) - @timeit SimMetaData.HourGlass "01 Calculate IndexCounter" begin - SimMetaData.Δx = UpdateΔx!(SimMetaData.Δx, Positionₙ⁺, SimParticles.Position) - ShouldRebuild = SimMetaData.Δx >= SimKernel.h - - if ShouldRebuild - @timeit SimMetaData.HourGlass "01a Actual Calculate IndexCounter" begin - SimMetaData.IndexCounter = UpdateNeighbors!(SimParticles, SimKernel.H⁻¹, ParticleRanges, UniqueCells, CellLookup, ParticleOrder, CellOffsets) - end - SimMetaData.Δx = zero(eltype(dρdtI)) - UniqueCellsView = view(UniqueCells, 1:SimMetaData.IndexCounter) - BuildNeighborCellLists!(NeighborCellLists, FullStencil, UniqueCellsView, ParticleRanges, CellLookup) - end + @timeit SimMetaData.HourGlass "01 Calculate IndexCounter" SimMetaData.Δx = UpdateΔxCUDA!(SimMetaData.Δx, CUDASupport, CUDAParticles, Launch) + if SimMetaData.Δx >= SimKernel.h + @timeit SimMetaData.HourGlass "01a Rebuild Neighbor Grid" RebuildNeighborGridCUDA!(CUDAGrid, CUDAParticles, SimKernel, Launch) + SimMetaData.Δx = zero(eltype(CUDASupport.DρdtI)) end - @timeit SimMetaData.HourGlass "Motion" ProgressMotion(SimParticles, dt₂, MotionDefinition, SimMetaData) - - @timeit SimMetaData.HourGlass "02 Pressure" Pressure!(SimParticles.Pressure, SimParticles.Density, SimConstants) - if SimMetaData isa SimulationMetaData{Dimensions, FloatType, SMode, KMode, NoMDBC, LMode} where {SMode, KMode, LMode} - @timeit SimMetaData.HourGlass "03 Apply MDBC before Half TimeStep" SPHCellList.ApplyMDBCBeforeHalf!(SimMetaData) - else - @timeit SimMetaData.HourGlass "03 Apply MDBC before Half TimeStep" SPHCellList.ApplyMDBCBeforeHalf!( - SimMetaData, SimKernel, SimConstants, SimParticles, ParticleRanges, CellLookup, Position, - Density, GhostPoints, GhostNormals, ParticleType; ParticleOrder = ParticleOrder, - ) + if MotionBuffers !== nothing + @timeit SimMetaData.HourGlass "Motion" ProgressMotionCUDA!(CUDAParticles, MotionBuffers, SimMetaData.TotalTime, dt₂, Launch) end + @timeit SimMetaData.HourGlass "02 Pressure" PressureCUDA!(CUDAParticles.Pressure, CUDAParticles.Density, SimConstants, Launch) + @timeit SimMetaData.HourGlass "04 First NeighborLoop" NeighborLoopPerParticleCUDA!( - SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData, - SimConstants, SimParticles, ParticleRanges, CellLookup, - NeighborCellLists, dρdtI, Acceleration, ∇Cᵢ, ∇◌rᵢ, - ParticleOrder = ParticleOrder, + SimDensityDiffusion, SimViscosity, SimKernel, SimConstants, + CUDAParticles, CUDAGrid, CUDASupport.DρdtI, CUDAParticles.Acceleration, Launch, ) - @timeit SimMetaData.HourGlass "05 Update To Half TimeStep" HalfTimeStep(SimMetaData, SimConstants, SimParticles, Positionₙ⁺, Velocityₙ⁺, ρₙ⁺, dρdtI, dt₂) + @timeit SimMetaData.HourGlass "05 Update To Half TimeStep" HalfTimeStepCUDA!(CUDAParticles, CUDASupport, dt₂, SimConstants, Launch) - @timeit SimMetaData.HourGlass "06 Half LimitDensityAtBoundary" LimitDensityAtBoundary!(ρₙ⁺, SimConstants.ρ₀, MotionLimiter) + @timeit SimMetaData.HourGlass "06 Half LimitDensityAtBoundary" LimitDensityAtBoundaryCUDA!(CUDASupport.ρₙ⁺, SimConstants.ρ₀, CUDAParticles.MotionLimiter, Launch) - @timeit SimMetaData.HourGlass "Motion" ProgressMotion(SimParticles, dt₂, MotionDefinition, SimMetaData) + if MotionBuffers !== nothing + @timeit SimMetaData.HourGlass "Motion" ProgressMotionCUDA!(CUDAParticles, MotionBuffers, SimMetaData.TotalTime, dt₂, Launch) + end - @timeit SimMetaData.HourGlass "07 Pressure" Pressure!(SimParticles.Pressure, ρₙ⁺, SimConstants) + @timeit SimMetaData.HourGlass "07 Pressure" PressureCUDA!(CUDAParticles.Pressure, CUDASupport.ρₙ⁺, SimConstants, Launch) @timeit SimMetaData.HourGlass "08 Second NeighborLoop" NeighborLoopPerParticleCUDA!( - SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData, - SimConstants, SimParticles, ParticleRanges, CellLookup, - NeighborCellLists, dρdtI, Acceleration, ∇Cᵢ, ∇◌rᵢ, - max_visc, min_dt_force, - ParticleOrder = ParticleOrder, - Position = Positionₙ⁺, - Density = ρₙ⁺, - Velocity = Velocityₙ⁺, + SimDensityDiffusion, SimViscosity, SimKernel, SimConstants, + CUDAParticles, CUDAGrid, CUDASupport.DρdtI, CUDAParticles.Acceleration, Launch, + Position = CUDASupport.Positionₙ⁺, + Density = CUDASupport.ρₙ⁺, + Velocity = CUDASupport.Velocityₙ⁺, ) - @timeit SimMetaData.HourGlass "09 Final LimitDensityAtBoundary" LimitDensityAtBoundary!(Density, SimConstants.ρ₀, MotionLimiter) + @timeit SimMetaData.HourGlass "09 Final LimitDensityAtBoundary" LimitDensityAtBoundaryCUDA!(CUDAParticles.Density, SimConstants.ρ₀, CUDAParticles.MotionLimiter, Launch) - @timeit SimMetaData.HourGlass "10 Final Density" DensityEpsi!(Density, dρdtI, ρₙ⁺, dt) + @timeit SimMetaData.HourGlass "10 Final Density" DensityEpsiCUDA!(CUDAParticles.Density, CUDASupport.DρdtI, CUDASupport.ρₙ⁺, dt, Launch) - @timeit SimMetaData.HourGlass "11 Update To Final TimeStep" FullTimeStep(SimMetaData, SimKernel, SimConstants, SimParticles, ∇Cᵢ, ∇◌rᵢ, dt) + @timeit SimMetaData.HourGlass "11 Update To Final TimeStep" FullTimeStepCUDA!(CUDAParticles, dt, SimConstants, Launch) @timeit SimMetaData.HourGlass "12 Update MetaData" UpdateMetaData!(SimMetaData, dt) - @timeit SimMetaData.HourGlass "13 Update TimeStep" dt = FinalizeTimeStep(max_visc, min_dt_force, SimConstants, SimKernel) + @timeit SimMetaData.HourGlass "13 Update TimeStep" dt = ΔtCUDA(CUDASupport, CUDAParticles, SimConstants, SimKernel, Launch) + dt₂ = dt * 0.5 end + SyncParticlesToHost!(SimParticles, CUDAParticles) + return nothing end @@ -380,8 +690,8 @@ function RunSimulationCUDA(;SimGeometry::Vector{Geometry{Dimensions, FloatType}} if !CUDAAvailable() error("CUDA is not available; cannot run RunSimulationCUDA.") end - if !(SimMetaData isa SimulationMetaData{Dimensions, FloatType, NoShifting, NoKernelOutput, BMode, LMode} where {BMode, LMode}) - error("RunSimulationCUDA currently supports NoShifting and NoKernelOutput configurations.") + if !(SimMetaData isa SimulationMetaData{Dimensions, FloatType, NoShifting, NoKernelOutput, NoMDBC, LMode} where {LMode}) + error("RunSimulationCUDA currently supports NoShifting, NoKernelOutput, and NoMDBC configurations.") end NumberOfPoints = length(SimParticles) @@ -423,14 +733,18 @@ function RunSimulationCUDA(;SimGeometry::Vector{Geometry{Dimensions, FloatType}} end MotionDefinition = SPHCellList.GenerateMotionDetails(SimParticles, SimGeometry, Dimensions, FloatType) + CUDAParticles = BuildCUDAParticleBuffers(SimParticles) + CUDASupport = BuildCUDASupportBuffers(dρdtI, Velocityₙ⁺, Positionₙ⁺, ρₙ⁺, ∇Cᵢ, ∇◌rᵢ) + MotionBuffers = MotionDefinition === nothing ? nothing : BuildCUDAMotionBuffers(SimParticles, MotionDefinition, Val(Dimensions), FloatType) + CUDAGrid = BuildCUDAGridBuffers(SimParticles, FullStencil, Val(Dimensions)) @inbounds while true @timeit SimMetaData.HourGlass "00 SimulationLoop" SimulationLoopCUDA( SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData, SimConstants, SimParticles, FullStencil, ParticleRanges, UniqueCells, CellLookup, ParticleOrder, CellOffsets, - NeighborCellLists, dρdtI, Velocityₙ⁺, Positionₙ⁺, ρₙ⁺, - ∇Cᵢ, ∇◌rᵢ, MotionDefinition, + NeighborCellLists, CUDAParticles, CUDASupport, + MotionBuffers, CUDAGrid, ) push!(SimMetaData.TimeSteps, SimMetaData.CurrentTimeStep) @@ -438,6 +752,13 @@ function RunSimulationCUDA(;SimGeometry::Vector{Geometry{Dimensions, FloatType}} SimMetaData.OutputIterationCounter += 1 + if SimMetaData.ExportGridCellParticleCounts + SimMetaData.IndexCounter = UpdateNeighbors!(SimParticles, SimKernel.H⁻¹, ParticleRanges, UniqueCells, CellLookup, ParticleOrder, CellOffsets) + if SimMetaData.IndexCounter > 0 + unique_cells_view = view(UniqueCells, 1:SimMetaData.IndexCounter) + BuildNeighborCellLists!(NeighborCellLists, FullStencil, unique_cells_view, ParticleRanges, CellLookup) + end + end UniqueCellsView = view(UniqueCells, 1:SimMetaData.IndexCounter) @timeit SimMetaData.HourGlass "13 Determine Output" SPHCellList.DetermineOutput( diff --git a/src/SPHCellList.jl b/src/SPHCellList.jl index 7d149e2e..d3b721bb 100644 --- a/src/SPHCellList.jl +++ b/src/SPHCellList.jl @@ -39,8 +39,7 @@ using LinearAlgebra SimConstants, SimParticles, ParticleRanges, CellLookup, NeighborCellLists, dρdtI, Acceleration, ∇Cᵢ, - ∇◌rᵢ, max_visc = nothing, - min_dt_force = nothing; + ∇◌rᵢ; Position = SimParticles.Position, Density = SimParticles.Density, Pressure = SimParticles.Pressure, @@ -84,7 +83,6 @@ using LinearAlgebra dρdtI[i] = dρdt_acc Acceleration[i] = acc_acc - UpdateTimeStepBuffers!(max_visc, min_dt_force, i, Position[i], Velocity[i], acc_acc, SimKernel) end return nothing @@ -95,8 +93,7 @@ using LinearAlgebra SimConstants, SimParticles, ParticleRanges, CellLookup, NeighborCellLists, dρdtI, Acceleration, ∇Cᵢ, - ∇◌rᵢ, max_visc = nothing, - min_dt_force = nothing; + ∇◌rᵢ; Position = SimParticles.Position, Density = SimParticles.Density, Pressure = SimParticles.Pressure, @@ -149,7 +146,6 @@ using LinearAlgebra Acceleration[i] = acc_acc Kernel[i] = kernel_acc KernelGradient[i] = kernel_grad_acc - UpdateTimeStepBuffers!(max_visc, min_dt_force, i, Position[i], Velocity[i], acc_acc, SimKernel) end return nothing @@ -160,8 +156,7 @@ using LinearAlgebra SimConstants, SimParticles, ParticleRanges, CellLookup, NeighborCellLists, dρdtI, Acceleration, ∇Cᵢ, - ∇◌rᵢ, max_visc = nothing, - min_dt_force = nothing; + ∇◌rᵢ; Position = SimParticles.Position, Density = SimParticles.Density, Pressure = SimParticles.Pressure, @@ -213,7 +208,6 @@ using LinearAlgebra Acceleration[i] = acc_acc ∇Cᵢ[i] = shift_c_acc ∇◌rᵢ[i] = shift_r_acc - UpdateTimeStepBuffers!(max_visc, min_dt_force, i, Position[i], Velocity[i], acc_acc, SimKernel) end return nothing @@ -224,8 +218,7 @@ using LinearAlgebra SimConstants, SimParticles, ParticleRanges, CellLookup, NeighborCellLists, dρdtI, Acceleration, ∇Cᵢ, - ∇◌rᵢ, max_visc = nothing, - min_dt_force = nothing; + ∇◌rᵢ; Position = SimParticles.Position, Density = SimParticles.Density, Pressure = SimParticles.Pressure, @@ -283,8 +276,6 @@ using LinearAlgebra KernelGradient[i] = kernel_grad_acc ∇Cᵢ[i] = shift_c_acc ∇◌rᵢ[i] = shift_r_acc - UpdateTimeStepBuffers!(max_visc, min_dt_force, i, Position[i], - Velocity[i], acc_acc, SimKernel) end return nothing @@ -729,9 +720,6 @@ using LinearAlgebra dt = Δt(Position, Velocity, Acceleration, SimConstants, SimKernel) @no_escape begin - max_visc = @alloc(FloatType, length(Position)) - min_dt_force = @alloc(FloatType, length(Position)) - dt₂ = dt * 0.5 while SimMetaData.TotalTime <= next_output_time(SimMetaData) @@ -795,7 +783,6 @@ using LinearAlgebra SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData, SimConstants, SimParticles, ParticleRanges, CellLookup, NeighborCellLists, dρdtI, Acceleration, ∇Cᵢ, ∇◌rᵢ, - max_visc, min_dt_force, ParticleOrder = ParticleOrder, Position = Positionₙ⁺, Density = ρₙ⁺, @@ -806,7 +793,6 @@ using LinearAlgebra SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData, SimConstants, SimParticles, ParticleRanges, CellLookup, NeighborCellLists, dρdtI, Acceleration, ∇Cᵢ, ∇◌rᵢ, - max_visc, min_dt_force, ParticleOrder = ParticleOrder, Position = Positionₙ⁺, Density = ρₙ⁺, @@ -822,7 +808,8 @@ using LinearAlgebra @timeit SimMetaData.HourGlass "12 Update MetaData" UpdateMetaData!(SimMetaData, dt) - @timeit SimMetaData.HourGlass "13 Update TimeStep" dt = FinalizeTimeStep(max_visc, min_dt_force, SimConstants, SimKernel) + @timeit SimMetaData.HourGlass "13 Update TimeStep" dt = Δt(Positionₙ⁺, Velocityₙ⁺, Acceleration, SimConstants, SimKernel) + dt₂ = dt * 0.5 end end diff --git a/src/TimeStepping.jl b/src/TimeStepping.jl index 0c795a95..e3b22a48 100644 --- a/src/TimeStepping.jl +++ b/src/TimeStepping.jl @@ -1,6 +1,6 @@ module TimeStepping -export Δt, FinalizeTimeStep, UpdateTimeStepBuffers!, next_output_time, ProgressMotion, HalfTimeStep, FullTimeStep +export Δt, next_output_time, ProgressMotion, HalfTimeStep, FullTimeStep using LinearAlgebra using Parameters @@ -10,42 +10,6 @@ using ..SimulationEquations using ..SimulationGeometry using ..SimulationMetaDataConfiguration -@inline function UpdateTimeStepBuffers!(::Nothing, ::Nothing, index, position, - velocity, acceleration, sim_kernel) - return nothing -end - -@inline function UpdateTimeStepBuffers!(max_visc, min_dt_force, index, position, - velocity, acceleration, sim_kernel) - h = sim_kernel.h - η² = sim_kernel.η² - r_sq = sqrt(dot(position, position))^2 - curr_visc = abs(h * dot(velocity, position) / (r_sq + η²)) - a_mag = norm(acceleration) - curr_dt_force = a_mag > 0 ? sqrt(h / a_mag) : typemax(eltype(min_dt_force)) - @inbounds begin - max_visc[index] = curr_visc - min_dt_force[index] = curr_dt_force - end - return nothing -end - -""" - FinalizeTimeStep(max_visc, min_dt_force, SimulationConstants, SPHKernel) - -Compute the CFL-limited time step from the per-particle buffers. -""" -function FinalizeTimeStep(max_visc, min_dt_force, SimulationConstants, SPHKernel) - @unpack c₀, CFL = SimulationConstants - @unpack h = SPHKernel - - global_visc = maximum(max_visc) - global_dt_force = minimum(min_dt_force) - - dt2 = h / (c₀ + global_visc) - return CFL * min(global_dt_force, dt2) -end - """ Δt(Position, Velocity, Acceleration, SimulationConstants, SPHKernel)