Skip to content

Make half-step wrapper views array-like and integrate into time-stepping #154

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 1 commit into
base: codex/refactor-to-avoid-if-statement
Choose a base branch
from
Open

Make half-step wrapper views array-like and integrate into time-stepping #154

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
11 changes: 4 additions & 7 deletions example/StillWedgeMDBC.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -7,9 +7,10 @@ let
SimConstantsWedge = SimulationConstants{FloatType}(dx=0.02,c₀=42.48576250492629, δᵩ = 0.1, CFL=0.5)
# SimConstantsWedge = SimulationConstants{FloatType}(dx=0.01,c₀=43.4, δᵩ = 0.1, CFL=0.2)
#
OutputDirectory = joinpath(pwd(), "output", "StillWedge2D_MDBC")
SimMetaDataWedge = SimulationMetaData{Dimensions,FloatType,NoShifting,NoKernelOutput,SimpleMDBC,StoreLog}(
SimulationName="StillWedge",
SaveLocation="W:/Simulations/StillWedge2D_MDBC",
SaveLocation=OutputDirectory,
SimulationTime=4.0,
OutputTimes=0.01,
VisualizeInParaview=true,
Expand All @@ -19,9 +20,7 @@ let
)

# If save directory is not already made, make it
if !isdir(SimMetaDataWedge.SaveLocation)
mkdir(SimMetaDataWedge.SaveLocation)
end
mkpath(SimMetaDataWedge.SaveLocation)

# Assuming SimConstantsWedge is defined somewhere else with the field `dx`
FixedBoundary = Geometry{Dimensions, FloatType}(
Expand Down Expand Up @@ -49,7 +48,7 @@ let

SimKernel = SPHKernelInstance{Dimensions, FloatType}(WendlandC2(); dx = SimConstantsWedge.dx)

@profview RunSimulation(
RunSimulation(
SimGeometry = SimulationGeometry,
SimMetaData = SimMetaDataWedge,
SimConstants = SimConstantsWedge,
Expand Down Expand Up @@ -97,5 +96,3 @@ let

# display(plt)
end


93 changes: 73 additions & 20 deletions src/SPHCellList.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -32,6 +32,58 @@ using Base.Threads
using LinearAlgebra
using Bumper

struct HalfStepPosition{P, V, M, T, PT} <: AbstractVector{PT}
Position::P
Velocity::V
MotionLimiter::M
Δt₂::T
end

HalfStepPosition(Position, Velocity, MotionLimiter, Δt₂) = HalfStepPosition(Position, Velocity, MotionLimiter, Δt₂, eltype(Position))

Base.IndexStyle(::Type{<:HalfStepPosition}) = IndexLinear()
Base.axes(half::HalfStepPosition) = axes(half.Position)
Base.size(half::HalfStepPosition) = size(half.Position)
Base.length(half::HalfStepPosition) = length(half.Position)
Base.getindex(half::HalfStepPosition, i::Int) = half.Position[i] + half.Velocity[i] * half.Δt₂ * half.MotionLimiter[i]

struct HalfStepVelocity{V, A, M, T, VT} <: AbstractVector{VT}
Velocity::V
Acceleration::A
MotionLimiter::M
Δt₂::T
end

HalfStepVelocity(Velocity, Acceleration, MotionLimiter, Δt₂) = HalfStepVelocity(Velocity, Acceleration, MotionLimiter, Δt₂, eltype(Velocity))

Base.IndexStyle(::Type{<:HalfStepVelocity}) = IndexLinear()
Base.axes(half::HalfStepVelocity) = axes(half.Velocity)
Base.size(half::HalfStepVelocity) = size(half.Velocity)
Base.length(half::HalfStepVelocity) = length(half.Velocity)
Base.getindex(half::HalfStepVelocity, i::Int) = half.Velocity[i] + half.Acceleration[i] * half.Δt₂ * half.MotionLimiter[i]

struct HalfStepDensity{D, R, M, T, DT} <: AbstractVector{DT}
Density::D
dρdtI::R
MotionLimiter::M
ρ₀::T
Δt₂::T
end

HalfStepDensity(Density, dρdtI, MotionLimiter, ρ₀, Δt₂) = HalfStepDensity(Density, dρdtI, MotionLimiter, ρ₀, Δt₂, eltype(Density))

Base.IndexStyle(::Type{<:HalfStepDensity}) = IndexLinear()
Base.axes(half::HalfStepDensity) = axes(half.Density)
Base.size(half::HalfStepDensity) = size(half.Density)
Base.length(half::HalfStepDensity) = length(half.Density)
function Base.getindex(half::HalfStepDensity, i::Int)
density = half.Density[i] + half.dρdtI[i] * half.Δt₂
if (density < half.ρ₀) * !Bool(half.MotionLimiter[i])
density = half.ρ₀
end
return density
end

function NeighborLoopPerParticle!(SimDensityDiffusion::SDD, SimViscosity::SV, SimKernel,
SimMetaData::SimulationMetaData{D,T,NoShifting,NoKernelOutput,B,L},
SimConstants, SimParticles, ParticleRanges,
Expand Down Expand Up @@ -688,9 +740,7 @@ using LinearAlgebra
SimMetaData::SimulationMetaData{Dimensions, FloatType, SMode, KMode, BMode, LMode},
SimConstants, SimParticles, FullStencil,
ParticleRanges, UniqueCells, CellDict,
SortingScratchSpace,
NeighborCellLists, dρdtI, Velocityₙ⁺,
Positionₙ⁺, ρₙ⁺, ∇Cᵢ, ∇◌rᵢ,
SortingScratchSpace, NeighborCellLists,
MotionDefinition::Union{
Nothing,
AbstractVector{
Expand All @@ -717,14 +767,25 @@ using LinearAlgebra
dt = SimConstants.CFL * (SimKernel.h / SimConstants.c₀)

@no_escape begin
AccelerationMax = @alloc(FloatType, length(Position))
PositionType = eltype(Position)
PositionUnderlyingType = eltype(PositionType)
AccelerationMax = similar(Density)
dρdtI = similar(Density)
∇Cᵢ = SMode <: NoShifting ? Vector{PositionType}(undef, 0) : similar(Position)
∇◌rᵢ = SMode <: NoShifting ? Vector{PositionUnderlyingType}(undef, 0) : similar(Density)
dt₂ = dt * 0.5
fill!(dρdtI, zero(FloatType))
if !(SMode <: NoShifting)
fill!(∇Cᵢ, zero(PositionType))
fill!(∇◌rᵢ, zero(PositionUnderlyingType))
end

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
PositionHalfStep = HalfStepPosition(Position, Velocity, MotionLimiter, dt₂)
SimMetaData.Δx = UpdateΔx!(SimMetaData.Δx, PositionHalfStep, SimParticles.Position)
ShouldRebuild = SimMetaData.Δx >= SimKernel.h || SimMetaData.IndexCounter == 0

# println("Δx: ", Δx, "h: ", SimKernel.h," dt: ", SimMetaData.CurrentTimeStep, " Iteration: ", SimMetaData.Iteration, " TotalTime: ", SimMetaData.TotalTime, " OutputIterationCounter: ", SimMetaData.OutputIterationCounter)

Expand Down Expand Up @@ -757,26 +818,21 @@ using LinearAlgebra
)


@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.ρ₀, MotionLimiter)

@timeit SimMetaData.HourGlass "Motion" ProgressMotion(SimParticles, dt₂, MotionDefinition, SimMetaData)

@timeit SimMetaData.HourGlass "07 Pressure" Pressure!(SimParticles.Pressure, ρₙ⁺,SimConstants)
@timeit SimMetaData.HourGlass "07 Pressure" PressureHalfStep!(SimParticles.Pressure, Density, dρdtI, MotionLimiter, SimConstants.ρ₀, dt₂, SimConstants)
@timeit SimMetaData.HourGlass "08 Second NeighborLoop" NeighborLoopPerParticle!(
SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData,
SimConstants, SimParticles, ParticleRanges, CellDict,
NeighborCellLists, dρdtI, Acceleration, ∇Cᵢ, ∇◌rᵢ, AccelerationMax,
Position = Positionₙ⁺,
Density = ρₙ⁺,
Velocity = Velocityₙ⁺,
Position = HalfStepPosition(Position, Velocity, MotionLimiter, dt₂),
Density = HalfStepDensity(Density, dρdtI, MotionLimiter, SimConstants.ρ₀, dt₂),
Velocity = HalfStepVelocity(Velocity, Acceleration, MotionLimiter, dt₂),
)

@timeit SimMetaData.HourGlass "09 Final LimitDensityAtBoundary" LimitDensityAtBoundary!(Density, SimConstants.ρ₀, MotionLimiter)

@timeit SimMetaData.HourGlass "10 Final Density" DensityEpsi!(Density, dρdtI, ρₙ⁺, dt)
@timeit SimMetaData.HourGlass "10 Final Density" DensityEpsi!(Density, dρdtI, MotionLimiter, SimConstants.ρ₀, dt, dt₂)

@timeit SimMetaData.HourGlass "11 Update To Final TimeStep" FullTimeStep(SimMetaData, SimKernel, SimConstants, SimParticles, ∇Cᵢ, ∇◌rᵢ, dt)

Expand Down Expand Up @@ -805,8 +861,6 @@ using LinearAlgebra
) where {Dimensions,FloatType,SMode,KMode,BMode,LMode,SV<:SPHViscosity,SDD<:SPHDensityDiffusion}

NumberOfPoints = length(SimParticles)

dρdtI, Velocityₙ⁺, Positionₙ⁺, ρₙ⁺, ∇Cᵢ, ∇◌rᵢ = AllocateSupportDataStructures(SimMetaData, SimParticles.Position)

LoadMDBCNormals!(SimMetaData, SimParticles, ParticleNormalsPath)

Expand Down Expand Up @@ -859,8 +913,7 @@ using LinearAlgebra
SimDensityDiffusion, SimViscosity, SimKernel, SimMetaData,
SimConstants, SimParticles, FullStencil, ParticleRanges,
UniqueCells, CellDict, SortingScratchSpace,
NeighborCellLists, dρdtI, Velocityₙ⁺, Positionₙ⁺, ρₙ⁺,
∇Cᵢ, ∇◌rᵢ, MotionDefinition,
NeighborCellLists, MotionDefinition,
)
push!(SimMetaData.TimeSteps, SimMetaData.CurrentTimeStep)

Expand Down
22 changes: 21 additions & 1 deletion src/SPHNeighborList.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3,6 +3,7 @@ module SPHNeighborList
export ConstructStencil, ExtractCells!, UpdateNeighbors!, BuildNeighborCellLists!, ComputeCellParticleCounts, ComputeCellNeighborCounts, UpdateΔx!

using StaticArrays
using LinearAlgebra: dot

function ConstructStencil(V::Val{d}) where d
return CartesianIndices(ntuple(_ -> -1:1, V))
Expand Down Expand Up @@ -77,6 +78,10 @@ Updates the neighbor list and sorts particles by their cell indices.
"""
function UpdateNeighbors!(Particles, InverseCutOff, SortingScratchSpace,
ParticleRanges, UniqueCells, CellDict)
RequiredLen = length(Particles) + 2
if length(ParticleRanges) < RequiredLen
resize!(ParticleRanges, RequiredLen)
end
ExtractCells!(Particles, InverseCutOff)

sort!(Particles, by = p -> p.Cells; scratch=SortingScratchSpace)
Expand All @@ -97,7 +102,7 @@ function UpdateNeighbors!(Particles, InverseCutOff, SortingScratchSpace,
CellDict[Cells[Index]] = IndexCounter
end
end
ParticleRanges[IndexCounter + 1] = length(ParticleRanges)
ParticleRanges[IndexCounter + 1] = length(Particles) + 1

return IndexCounter
end
Expand Down Expand Up @@ -149,4 +154,19 @@ Returns the new Δx.
return Δx + 4 * maxd
end

@inline function UpdateΔx!(Δx::T,
posₙ⁺,
pos::AbstractVector{SVector{D, T}}) where {D, T<:Real}
maxd = zero(T)
@inbounds for i in eachindex(pos)
diff = posₙ⁺[i] - pos[i]
sumsq = dot(diff, diff)
nrm = sqrt(sumsq)
if nrm > maxd
maxd = nrm
end
end
return Δx + 4 * maxd
end

end
24 changes: 23 additions & 1 deletion src/SimulationEquations.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
module SimulationEquations

export EquationOfState, EquationOfStateGamma7, Pressure!, DensityEpsi!, LimitDensityAtBoundary!, ConstructGravitySVector, InverseHydrostaticEquationOfState, Estimate7thRoot
export EquationOfState, EquationOfStateGamma7, Pressure!, PressureHalfStep!, DensityEpsi!, LimitDensityAtBoundary!, ConstructGravitySVector, InverseHydrostaticEquationOfState, Estimate7thRoot

using StaticArrays
using Parameters
Expand All @@ -23,6 +23,17 @@ end
end
end

@inline function PressureHalfStep!(Press, Density, dρdtI, MotionLimiter, ρ₀, Δt₂, SimulationConstants)
@unpack c₀ = SimulationConstants
@inbounds for i ∈ eachindex(Press, Density, dρdtI, MotionLimiter)
ρₙ⁺ = Density[i] + dρdtI[i] * Δt₂
if (ρₙ⁺ < ρ₀) * !Bool(MotionLimiter[i])
ρₙ⁺ = ρ₀
end
Press[i] = EquationOfStateGamma7(ρₙ⁺, c₀, ρ₀)
end
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)
Expand All @@ -32,6 +43,17 @@ end
end
end

@inline function DensityEpsi!(Density, dρdtI, MotionLimiter, ρ₀, Δt, Δt₂)
@inbounds for i in eachindex(Density, dρdtI, MotionLimiter)
ρₙ⁺ = Density[i] + dρdtI[i] * Δt₂
if (ρₙ⁺ < ρ₀) * !Bool(MotionLimiter[i])
ρₙ⁺ = ρ₀
end
epsi = - (dρdtI[i] / ρₙ⁺) * Δt
Density[i] *= (2 - epsi) / (2 + epsi)
end
end

# This version of the function using !Bool(MotionLimiter) instead of BoundaryBool
@inline function LimitDensityAtBoundary!(Density,ρ₀, MotionLimiter)
@inbounds for i in eachindex(Density)
Expand Down