From 980b021454d53b67ddc22cb1840db9e222597829 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Thu, 21 Nov 2024 15:35:38 +0100 Subject: [PATCH] fix splitting --- src/refinement/refinement.jl | 29 ++++++++++++++-------------- src/refinement/refinement_pattern.jl | 21 ++++++-------------- src/refinement/split.jl | 28 +++++++++++++-------------- 3 files changed, 35 insertions(+), 43 deletions(-) diff --git a/src/refinement/refinement.jl b/src/refinement/refinement.jl index 616d3e98d..58beef74f 100644 --- a/src/refinement/refinement.jl +++ b/src/refinement/refinement.jl @@ -4,15 +4,14 @@ include("split.jl") include("merge.jl") struct ParticleRefinement{SP, RC, ELTYPE} - refinement_pattern :: SP - refinement_criteria :: RC - max_spacing_ratio :: ELTYPE - mass_ref :: Vector{ELTYPE} # length(mass_ref) == nparticles - merge_candidates :: Vector{Int} # length(merge_candidates) == nparticles - delete_candidates :: Vector{Bool} # length(delete_candidates) == nparticles - split_candidates :: Vector{Int} # variable length - n_particles_before_resize :: Ref{Int} - n_new_particles :: Ref{Int} + refinement_pattern :: SP + refinement_criteria :: RC + max_spacing_ratio :: ELTYPE + mass_ref :: Vector{ELTYPE} # length(mass_ref) == nparticles + merge_candidates :: Vector{Int} # length(merge_candidates) == nparticles + delete_candidates :: Vector{Bool} # length(delete_candidates) == nparticles + split_candidates :: Vector{Int} # variable length + n_new_particles :: Ref{Int} end function ParticleRefinement(; refinement_pattern, max_spacing_ratio, @@ -25,7 +24,7 @@ function ParticleRefinement(; refinement_pattern, max_spacing_ratio, end return ParticleRefinement(refinement_pattern, refinement_criteria, max_spacing_ratio, - mass_ref, Int[], delete_candidates, Int[], Ref(0), Ref(0)) + mass_ref, Int[], delete_candidates, Int[], Ref(0)) end resize_refinement!(system) = system @@ -41,6 +40,8 @@ function resize_refinement!(system, particle_refinement) resize!(particle_refinement.delete_candidates, nparticles(system)) resize!(particle_refinement.merge_candidates, nparticles(system)) + particle_refinement.mass_ref .= system.mass + return system end @@ -160,9 +161,9 @@ function upate_smoothing_lengths!(system::FluidSystem, refinement, semi, u) end @inline function min_max_avg_spacing(system, semi, u_ode, system_coords, particle) - dp_min = Inf - dp_max = zero(eltype(system)) - dp_avg = zero(eltype(system)) + dp_min = particle_spacing(system, particle) + dp_max = particle_spacing(system, particle) + dp_avg = particle_spacing(system, particle) counter_neighbors = 0 foreach_system(semi) do neighbor_system @@ -183,7 +184,7 @@ end end end - dp_avg / counter_neighbors + dp_avg = counter_neighbors == 0 ? dp_avg : dp_avg / counter_neighbors return dp_min, dp_max, dp_avg end diff --git a/src/refinement/refinement_pattern.jl b/src/refinement/refinement_pattern.jl index 759fe8262..c01022314 100644 --- a/src/refinement/refinement_pattern.jl +++ b/src/refinement/refinement_pattern.jl @@ -3,6 +3,7 @@ struct CubicSplitting{ELTYPE} alpha :: ELTYPE center_particle :: Bool relative_position :: Vector{SVector{2, ELTYPE}} + n_children :: Int function CubicSplitting(; epsilon=0.5, alpha=0.5, center_particle=true) ELTYPE = typeof(epsilon) @@ -21,19 +22,16 @@ struct CubicSplitting{ELTYPE} push!(relative_position, SVector(zero(ELTYPE), zero(ELTYPE))) end - new{ELTYPE}(epsilon, alpha, center_particle, relative_position) + new{ELTYPE}(epsilon, alpha, center_particle, relative_position, 4 + center_particle) end end -@inline function nchilds(system, refinement_pattern::CubicSplitting) - return 4 + refinement_pattern.center_particle -end - struct TriangularSplitting{ELTYPE} epsilon :: ELTYPE alpha :: ELTYPE center_particle :: Bool relative_position :: Vector{SVector{2, ELTYPE}} + n_children :: Int function TriangularSplitting(; epsilon=0.5, alpha=0.5, center_particle=true) ELTYPE = typeof(epsilon) @@ -49,19 +47,16 @@ struct TriangularSplitting{ELTYPE} push!(relative_position, SVector(zero(ELTYPE), zero(ELTYPE))) end - new{ELTYPE}(epsilon, alpha, center_particle, relative_position) + new{ELTYPE}(epsilon, alpha, center_particle, relative_position, 3 + center_particle) end end -@inline function nchilds(system::System{2}, refinement_pattern::TriangularSplitting) - return 3 + refinement_pattern.center_particle -end - struct HexagonalSplitting{ELTYPE} epsilon :: ELTYPE alpha :: ELTYPE center_particle :: Bool relative_position :: Vector{SVector{2, ELTYPE}} + n_children :: Int function HexagonalSplitting(; epsilon=0.4, alpha=0.9, center_particle=true) ELTYPE = typeof(epsilon) @@ -84,10 +79,6 @@ struct HexagonalSplitting{ELTYPE} push!(relative_position, SVector(zero(ELTYPE), zero(ELTYPE))) end - new{ELTYPE}(epsilon, alpha, center_particle, relative_position) + new{ELTYPE}(epsilon, alpha, center_particle, relative_position, 6 + center_particle) end end - -@inline function nchilds(system::System{2}, refinement_pattern::HexagonalSplitting) - return 6 + refinement_pattern.center_particle -end diff --git a/src/refinement/split.jl b/src/refinement/split.jl index 47c83fef1..b7743d7d9 100644 --- a/src/refinement/split.jl +++ b/src/refinement/split.jl @@ -26,6 +26,7 @@ end @inline function collect_split_candidates!(system::FluidSystem, particle_refinement, v, u) (; mass_ref, max_spacing_ratio, split_candidates, refinement_pattern) = particle_refinement + (; n_children, center_particle) = refinement_pattern resize!(split_candidates, 0) @@ -38,10 +39,8 @@ end end end - n_childs_exclude_center = nchilds(system, refinement_pattern) - 1 - particle_refinement.n_new_particles[] = length(split_candidates) * - n_childs_exclude_center + (n_children - center_particle) return system end @@ -56,14 +55,15 @@ end @inline function split_particles!(system::FluidSystem, particle_refinement, v, u) (; smoothing_length) = system.cache - (; split_candidates, refinement_pattern, n_particles_before_resize) = particle_refinement - (; alpha, relative_position) = refinement_pattern + (; split_candidates, refinement_pattern, n_new_particles) = particle_refinement + (; alpha, relative_position, n_children, center_particle) = refinement_pattern + child_id_global = nparticles(system) - n_new_particles[] for particle in split_candidates smoothing_length_old = smoothing_length[particle] mass_old = system.mass[particle] - system.mass[particle] = mass_old / nchilds(system, refinement_pattern) + system.mass[particle] = mass_old / n_children p_a = particle_pressure(v, system, particle) rho_a = particle_density(v, system, particle) @@ -73,23 +73,23 @@ end pos_center = current_coords(u, system, particle) vel_center = current_velocity(v, system, particle) - for child_id_local in 1:(nchilds(system, refinement_pattern) - 1) - child = n_particles_before_resize[] + child_id_local + for child_id_local in 1:(n_children - center_particle) + child_id_global += 1 - system.mass[child] = mass_old / nchilds(system, refinement_pattern) + system.mass[child_id_global] = mass_old / n_children - set_particle_pressure!(v, system, child, p_a) + set_particle_pressure!(v, system, child_id_global, p_a) - set_particle_density!(v, system, child, rho_a) + set_particle_density!(v, system, child_id_global, rho_a) - smoothing_length[child] = alpha * smoothing_length_old + smoothing_length[child_id_global] = alpha * smoothing_length_old rel_pos = smoothing_length_old * relative_position[child_id_local] new_pos = pos_center + rel_pos for dim in 1:ndims(system) - u[dim, child] = new_pos[dim] - v[dim, child] = vel_center[dim] + u[dim, child_id_global] = new_pos[dim] + v[dim, child_id_global] = vel_center[dim] end end end