fix splitting

trixi-framework · Nov 21, 2024 · 980b021 · 980b021
1 parent e9da8f5
commit 980b021
Show file tree

Hide file tree

Showing 3 changed files with 35 additions and 43 deletions.
diff --git 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
@@ -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
@@ -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