From 61e147f53e4e4b7daae6c088115f9b2c0534b1f0 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Fri, 11 Apr 2025 17:24:27 +0200 Subject: [PATCH 1/3] add shifting for buffer particles --- src/callbacks/particle_shifting.jl | 93 +++++++++++++++++++- src/general/system.jl | 4 +- src/schemes/boundary/open_boundary/system.jl | 4 + 3 files changed, 97 insertions(+), 4 deletions(-) diff --git a/src/callbacks/particle_shifting.jl b/src/callbacks/particle_shifting.jl index 0a57f164c4..35122b2225 100644 --- a/src/callbacks/particle_shifting.jl +++ b/src/callbacks/particle_shifting.jl @@ -49,8 +49,7 @@ function particle_shifting!(u, v, system, v_ode, u_ode, semi, u_cache, dt) return u end -function particle_shifting!(u, v, system::FluidSystem, v_ode, u_ode, semi, - u_cache, dt) +function particle_shifting!(u, v, system::FluidSystem, v_ode, u_ode, semi, u_cache, dt) # Wrap the cache vector to an NDIMS x NPARTICLES matrix. # We need this buffer because we cannot safely update `u` while iterating over it. delta_r = wrap_u(u_cache, system, semi) @@ -121,3 +120,93 @@ function Base.show(io::IO, ::MIME"text/plain", summary_box(io, "ParticleShiftingCallback") end end + +# For `OpenBoundarySPHSystem` with an `InFlow`-zone, we need a high quality transition +# from the boundary zone to the fluid domain, as small perturbations can lead to instabilities. +function particle_shifting!(u, v, + system::OpenBoundarySPHSystem{<:Any, <:BoundaryZone{<:InFlow}}, + v_ode, u_ode, semi, u_cache, dt) + (; fluid_system, boundary_zone) = system + (; zone_origin, spanning_set) = boundary_zone + + # Wrap the cache vector to an NDIMS x NPARTICLES matrix. + # We need this buffer because we cannot safely update `u` while iterating over it. + delta_r = wrap_u(u_cache, system, semi) + set_zero!(delta_r) + + v_max = maximum(particle -> norm(current_velocity(v, system, particle)), + eachparticle(system)) + + # TODO this needs to be adapted to multi-resolution. + # Section 3.2 explains what else needs to be changed. + Wdx = smoothing_kernel(system, particle_spacing(system, 1), 1) + h = smoothing_length(system, 1) + + # We use the neighborhood search (NHS) of the fluid system. + # This is appropriate because we only want to shift the buffer particles near the transition plane. + # Therefore, we use half of the compact support as the maximum distance within which the buffer particles are shifted. + max_dist = 0.5 * compact_support(system_smoothing_kernel(system), + smoothing_length(system, 1)) + + foreach_system(semi) do neighbor_system + u_neighbor = wrap_u(u_ode, neighbor_system, semi) + v_neighbor = wrap_v(v_ode, neighbor_system, semi) + + neighborhood_search = get_neighborhood_search(fluid_system, neighbor_system, semi) + + @threaded semi for particle in each_moving_particle(system) + particle_coords = current_coords(u, system, particle) + + particle_position = particle_coords - zone_origin + + # Distance to transistion plane + dist = dot(particle_position, spanning_set[1]) + + if dist <= max_dist + for neighbor in PointNeighbors.eachneighbor(particle_coords, + neighborhood_search) + neighbor_coords = current_coords(u_neighbor, neighbor_system, neighbor) + + pos_diff = particle_coords - neighbor_coords + distance2 = dot(pos_diff, pos_diff) + + # Check if the neighbor is within the search radius + if distance2 <= neighborhood_search.search_radius^2 + distance = sqrt(distance2) + + m_b = hydrodynamic_mass(neighbor_system, neighbor) + rho_a = particle_density(v, system, particle) + rho_b = particle_density(v_neighbor, neighbor_system, neighbor) + + kernel = smoothing_kernel(system, distance, particle) + grad_kernel = smoothing_kernel_grad(system, pos_diff, distance, + particle) + + # According to p. 29 below Eq. 9 + R = 0.2 + n = 4 + + # Eq. 7 in Sun et al. (2017). + # CFL * Ma can be rewritten as Δt * v_max / h (see p. 29, right above Eq. 9). + delta_r_ = -dt * v_max * 4 * h * (1 + R * (kernel / Wdx)^n) * + m_b / (rho_a + rho_b) * grad_kernel + + # Write into the buffer + for i in eachindex(delta_r_) + @inbounds delta_r[i, particle] += delta_r_[i] + end + end + end + end + end + end + + # Add δ_r from the buffer to the current coordinates + @threaded semi for particle in eachparticle(system) + for i in axes(delta_r, 1) + @inbounds u[i, particle] += delta_r[i, particle] + end + end + + return u +end diff --git a/src/general/system.jl b/src/general/system.jl index e578f698b1..cb6b72fcee 100644 --- a/src/general/system.jl +++ b/src/general/system.jl @@ -103,8 +103,8 @@ end @inline set_particle_pressure!(v, system, particle, pressure) = v @inline function smoothing_kernel(system, distance, particle) - (; smoothing_kernel) = system - return kernel(smoothing_kernel, distance, smoothing_length(system, particle)) + return kernel(system_smoothing_kernel(system), distance, + smoothing_length(system, particle)) end @inline function smoothing_kernel_grad(system, pos_diff, distance, particle) diff --git a/src/schemes/boundary/open_boundary/system.jl b/src/schemes/boundary/open_boundary/system.jl index 9b2c3dc17b..115da0d62d 100644 --- a/src/schemes/boundary/open_boundary/system.jl +++ b/src/schemes/boundary/open_boundary/system.jl @@ -204,6 +204,10 @@ end update_callback_used!(system::OpenBoundarySPHSystem) = system.update_callback_used[] = true +function system_smoothing_kernel(system::OpenBoundarySPHSystem) + return system.fluid_system.smoothing_kernel +end + function smoothing_length(system::OpenBoundarySPHSystem, particle) return smoothing_length(system.fluid_system, particle) end From 9ee18bf800b25b7774ccd0df312339c999ce9aac Mon Sep 17 00:00:00 2001 From: LasNikas Date: Fri, 11 Apr 2025 17:56:23 +0200 Subject: [PATCH 2/3] typo --- src/callbacks/particle_shifting.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/callbacks/particle_shifting.jl b/src/callbacks/particle_shifting.jl index f946574919..71cc03789c 100644 --- a/src/callbacks/particle_shifting.jl +++ b/src/callbacks/particle_shifting.jl @@ -161,7 +161,7 @@ function particle_shifting!(u, v, particle_position = particle_coords - zone_origin - # Distance to transistion plane + # Distance to transition plane dist = dot(particle_position, spanning_set[1]) if dist <= max_dist From 683b09cdca579c775696def9932f4c429813373b Mon Sep 17 00:00:00 2001 From: LasNikas Date: Sun, 8 Jun 2025 23:35:54 +0200 Subject: [PATCH 3/3] rewrite --- src/callbacks/particle_shifting.jl | 170 +++++++++---------- src/schemes/boundary/open_boundary/system.jl | 4 + 2 files changed, 87 insertions(+), 87 deletions(-) diff --git a/src/callbacks/particle_shifting.jl b/src/callbacks/particle_shifting.jl index a74c852916..b6d175cc07 100644 --- a/src/callbacks/particle_shifting.jl +++ b/src/callbacks/particle_shifting.jl @@ -84,26 +84,8 @@ function particle_shifting!(u, v, system::FluidSystem, v_ode, u_ode, semi, semi; points=each_moving_particle(system)) do particle, neighbor, pos_diff, distance - m_b = hydrodynamic_mass(neighbor_system, neighbor) - rho_a = current_density(v, system, particle) - rho_b = current_density(v_neighbor, neighbor_system, neighbor) - - kernel = smoothing_kernel(system, distance, particle) - grad_kernel = smoothing_kernel_grad(system, pos_diff, distance, particle) - - # According to p. 29 below Eq. 9 - R = 2 // 10 - n = 4 - - # Eq. 7 in Sun et al. (2017). - # CFL * Ma can be rewritten as Δt * v_max / h (see p. 29, right above Eq. 9). - delta_r_ = -dt * v_max * 4 * h * (1 + R * (kernel / Wdx)^n) * - m_b / (rho_a + rho_b) * grad_kernel - - # Write into the buffer - for i in eachindex(delta_r_) - @inbounds delta_r[i, particle] += delta_r_[i] - end + apply_shifting!(delta_r, system, neighbor_system, v, v_neighbor, u, u_neighbor, + particle, neighbor, pos_diff, distance, v_max, Wdx, h, dt) end end @@ -117,48 +99,46 @@ function particle_shifting!(u, v, system::FluidSystem, v_ode, u_ode, semi, return u end -function Base.show(io::IO, cb::DiscreteCallback{<:Any, typeof(particle_shifting!)}) - @nospecialize cb # reduce precompilation time - print(io, "ParticleShiftingCallback()") +function particle_shifting!(u, v, system::OpenBoundarySPHSystem, v_ode, u_ode, semi, + vu_cache, dt) + particle_shifting!(u, v, system, system.boundary_zone, v_ode, u_ode, semi, + vu_cache, dt) end -function Base.show(io::IO, ::MIME"text/plain", - cb::DiscreteCallback{<:Any, typeof(particle_shifting!)}) - @nospecialize cb # reduce precompilation time - - if get(io, :compact, false) - show(io, cb) - else - summary_box(io, "ParticleShiftingCallback") - end +function particle_shifting!(u, v, system::OpenBoundarySPHSystem, ::BoundaryZone{OutFlow}, + v_ode, u_ode, semi, vu_cache, dt) + return system end -# For `OpenBoundarySPHSystem` with an `InFlow`-zone, we need a high quality transition -# from the boundary zone to the fluid domain, as small perturbations can lead to instabilities. -function particle_shifting!(u, v, - system::OpenBoundarySPHSystem{<:Any, <:BoundaryZone{<:InFlow}}, - v_ode, u_ode, semi, u_cache, dt) - (; fluid_system, boundary_zone) = system - (; zone_origin, spanning_set) = boundary_zone +function particle_shifting!(u, v, system::OpenBoundarySPHSystem, ::BoundaryZone{InFlow}, + v_ode, u_ode, semi, vu_cache, dt) + (; zone_origin, spanning_set) = system.boundary_zone # Wrap the cache vector to an NDIMS x NPARTICLES matrix. # We need this buffer because we cannot safely update `u` while iterating over it. + _, u_cache = vu_cache.x delta_r = wrap_u(u_cache, system, semi) set_zero!(delta_r) - v_max = maximum(particle -> norm(current_velocity(v, system, particle)), - eachparticle(system)) + # This has similar performance to `maximum(..., eachparticle(system))`, + # but is GPU-compatible. + v_max = maximum(x -> sqrt(dot(x, x)), + reinterpret(reshape, SVector{ndims(system), eltype(v)}, + current_velocity(v, system))) # TODO this needs to be adapted to multi-resolution. # Section 3.2 explains what else needs to be changed. Wdx = smoothing_kernel(system, particle_spacing(system, 1), 1) h = smoothing_length(system, 1) - # We use the neighborhood search (NHS) of the fluid system. - # This is appropriate because we only want to shift the buffer particles near the transition plane. - # Therefore, we use half of the compact support as the maximum distance within which the buffer particles are shifted. - max_dist = 0.5 * compact_support(system_smoothing_kernel(system), - smoothing_length(system, 1)) + fluid_system = corresponding_fluid_system(system, semi) + + max_dist2 = (compact_support(system_smoothing_kernel(system.fluid_system), + smoothing_length(system, 1)) / 2)^2 + + shift_candidates = findall(x -> dot(x - zone_origin, spanning_set[1]) <= max_dist2, + reinterpret(reshape, SVector{ndims(system), eltype(u)}, + active_coordinates(u, system))) foreach_system(semi) do neighbor_system u_neighbor = wrap_u(u_ode, neighbor_system, semi) @@ -166,49 +146,21 @@ function particle_shifting!(u, v, neighborhood_search = get_neighborhood_search(fluid_system, neighbor_system, semi) - @threaded semi for particle in each_moving_particle(system) - particle_coords = current_coords(u, system, particle) - - particle_position = particle_coords - zone_origin - - # Distance to transition plane - dist = dot(particle_position, spanning_set[1]) - - if dist <= max_dist - for neighbor in PointNeighbors.eachneighbor(particle_coords, - neighborhood_search) - neighbor_coords = current_coords(u_neighbor, neighbor_system, neighbor) - - pos_diff = particle_coords - neighbor_coords - distance2 = dot(pos_diff, pos_diff) - - # Check if the neighbor is within the search radius - if distance2 <= neighborhood_search.search_radius^2 - distance = sqrt(distance2) - - m_b = hydrodynamic_mass(neighbor_system, neighbor) - rho_a = particle_density(v, system, particle) - rho_b = particle_density(v_neighbor, neighbor_system, neighbor) - - kernel = smoothing_kernel(system, distance, particle) - grad_kernel = smoothing_kernel_grad(system, pos_diff, distance, - particle) - - # According to p. 29 below Eq. 9 - R = 0.2 - n = 4 + neighbor_coords = current_coordinates(u_neighbor, neighbor_system) - # Eq. 7 in Sun et al. (2017). - # CFL * Ma can be rewritten as Δt * v_max / h (see p. 29, right above Eq. 9). - delta_r_ = -dt * v_max * 4 * h * (1 + R * (kernel / Wdx)^n) * - m_b / (rho_a + rho_b) * grad_kernel + @threaded semi for particle in shift_candidates + particle_coords = current_coords(u, system, particle) - # Write into the buffer - for i in eachindex(delta_r_) - @inbounds delta_r[i, particle] += delta_r_[i] - end - end - end + # TODO: Not public API + PointNeighbors.foreach_neighbor(neighbor_coords, neighborhood_search, + particle, particle_coords, + neighborhood_search.search_radius) do particle, + neighbor, + pos_diff, + distance + apply_shifting!(delta_r, system, neighbor_system, v, v_neighbor, u, + u_neighbor, particle, neighbor, pos_diff, distance, + v_max, Wdx, h, dt) end end end @@ -222,3 +174,47 @@ function particle_shifting!(u, v, return u end + +function apply_shifting!(delta_r, system, neighbor_system, + v, v_neighbor, u, u_neighbor, + particle, neighbor, pos_diff, distance, v_max, Wdx, h, dt) + # Calculate the hydrodynamic mass and density + m_b = hydrodynamic_mass(neighbor_system, neighbor) + rho_a = current_density(v, system, particle) + rho_b = current_density(v_neighbor, neighbor_system, neighbor) + + kernel = smoothing_kernel(system, distance, particle) + grad_kernel = smoothing_kernel_grad(system, pos_diff, distance, particle) + + # According to p. 29 below Eq. 9 + R = 2 // 10 + n = 4 + + # Eq. 7 in Sun et al. (2017). + # CFL * Ma can be rewritten as Δt * v_max / h (see p. 29, right above Eq. 9). + delta_r_ = -dt * v_max * 4 * h * (1 + R * (kernel / Wdx)^n) * + m_b / (rho_a + rho_b) * grad_kernel + + # Write into the buffer + for i in eachindex(delta_r_) + @inbounds delta_r[i, particle] += delta_r_[i] + end + + return delta_r +end + +function Base.show(io::IO, cb::DiscreteCallback{<:Any, typeof(particle_shifting!)}) + @nospecialize cb # reduce precompilation time + print(io, "ParticleShiftingCallback()") +end + +function Base.show(io::IO, ::MIME"text/plain", + cb::DiscreteCallback{<:Any, typeof(particle_shifting!)}) + @nospecialize cb # reduce precompilation time + + if get(io, :compact, false) + show(io, cb) + else + summary_box(io, "ParticleShiftingCallback") + end +end diff --git a/src/schemes/boundary/open_boundary/system.jl b/src/schemes/boundary/open_boundary/system.jl index f299fd9286..8cefe5ef5b 100644 --- a/src/schemes/boundary/open_boundary/system.jl +++ b/src/schemes/boundary/open_boundary/system.jl @@ -236,6 +236,10 @@ function smoothing_length(system::OpenBoundarySPHSystem, particle) return system.smoothing_length end +function system_smoothing_kernel(system::OpenBoundarySPHSystem) + return system.fluid_system.smoothing_kernel +end + @inline hydrodynamic_mass(system::OpenBoundarySPHSystem, particle) = system.mass[particle] @inline function current_density(v, system::OpenBoundarySPHSystem)