add refinement criteria

src/multi_resolution/particle_refinement.jl

@@ -1,24 +1,10 @@
 struct ParticleRefinement{SP, ELTYPE}
     splitting_pattern       :: SP
-    particle_spacing_min    :: ELTYPE
-    smoothing_length_factor :: Vector{ELTYPE}
-end
-
-function initialize!(particle_refinement, semi)
-    (; particle_spacing_min) = particle_refinement
-
-    particle_spacing_min = zero(eltype(particle_refinement))
-
-    foreach_system(semi) do system
-        (; particle_spacing) = system.inital_condition
-        particle_spacing_min[] = min(particle_spacing_min, particle_spacing)
-    end
-
-    return particle_refinement
+    max_spacing_ratio       :: ELTYPE
 end
 
 function refinement!(semi, v_ode, u_ode, v_tmp, u_tmp, t)
-    # check refnement criteria
+    check_refinement_criteria!(semi, v_ode, u_ode)
 
     # Update the spacing of particles (Algorthm 1)
 
@@ -48,5 +34,87 @@ function refinement!(semi, v_ode, u_ode, v_tmp, u_tmp, t)
     return semi
 end
 
-function update_particle_spacing()
+function check_refinement_criteria!(semi::Semidiscretization, v_ode, u_ode)
+    foreach_system(semi) do system
+        check_refinement_criteria!(system, v_ode, u_ode)
+    end
+end
+
+@inline check_refinement_criteria!(system, v_ode, u_ode) = system
+
+@inline function check_refinement_criteria!(system::FluidSystem, v_ode, u_ode)
+    (; refinement_criteria) = system.particle_refinement
+    for criterion in refinement_criteria
+        criterion(system, semi, v, u)
+    end
+end
+
+function update_particle_spacing(semi::Semidiscretization, u_ode)
+    foreach_system(semi) do particle_system
+        u_particle_system = wrap_u(u_ode, particle_system, semi)
+        update_particle_spacing(particle_system, u_particle_system, semi)
+    end
+end
+
+@inline update_particle_spacing(system, u, semi) = system
+
+@inline function update_particle_spacing(system::FluidSystem, u, semi)
+    (; smoothing_length, smoothing_length_factor) = system.cache
+    system_coords = current_coordinates(u, system)
+
+    for particle in eachparticle(system)
+        dp_min, dp_max, dp_avg = min_max_avg_spacing(system, semi, system_coords, particle)
+
+        if dp_max / dp_min < max_spacing_ratio^3
+            new_spacing = min(dp_max, max_spacing_ratio * dp_min)
+        else
+            new_spacing = dp_avg
+        end
+
+        smoothing_length[particle] = smoothing_length_factor * new_spacing
+        system.mass[particle] = system.density[particle] * new_spacing^(ndims(system))
+    end
+
+    return system
+end
+
+@inline function min_max_avg_spacing(system, semi, system_coords, particle)
+    dp_min = Inf
+    dp_max = zero(eltype(system))
+    dp_avg = zero(eltype(system))
+    counter_neighbors = 0
+
+    foreach_system(semi) do neighbor_system
+        neighborhood_search = get_neighborhood_search(particle_system, neighbor_system,
+                                                      semi)
+
+        u_neighbor_system = wrap_u(u_ode, neighbor, semi)
+        neighbor_coords = current_coordinates(u_neighbor_system, neighbor_system)
+
+        PointNeighbors.foreach_neighbor(system_coords, neighbor_coords, neighborhood_search,
+                                        particle) do particle, neighbor, pos_diff, distance
+            dp_neighbor = particle_spacing(neighbor_system, neighbor)
+
+            dp_min = min(dp_min, dp_neighbor)
+            dp_max = max(dp_max, dp_neighbor)
+            dp_avg += dp_neighbor
+
+            counter_neighbors += 1
+        end
+    end
+
+    dp_avg / counter_neighbors
+
+    return dp_min, dp_max, dp_avg
+end
+
+@inline function particle_spacing(system, particle)
+    return particle_spacing(system, system.particle_refinement, particle)
+end
+
+@inline particle_spacing(system, ::Nohing, _) = system.initial_condition.particle_spacing
+
+@inline function particle_spacing(system, refinement, particle)
+    (; smoothing_length_factor) = system.cache
+    return smoothing_length(system, particle) / smoothing_length_factor
 end

src/multi_resolution/refinement_criteria.jl

@@ -0,0 +1,43 @@
+abstract type RefinementCriteria end
+
+struct SpatialRefinementCriterion <: RefinementCriteria end
+
+struct SolutionRefinementCriterion <: RefinementCriteria end
+
+@inline function (criterion::SpatialRefinementCriterion)(system, semi, v_ode, u_ode)
+    u = wrap_u(u_ode, system, semi)
+    system_coords = current_coordinates(u, system)
+
+    foreach_system(semi) do neighbor_system
+        set_particle_spacing!(system, neighbor_system, system_coords)
+    end
+    return system
+end
+
+@inline set_particle_spacing!(system, neighbor_system, system_coords) = system
+
+@inline function set_particle_spacing!(particle_system,
+                                       neighbor_system::Union{BoundarySPHSystem,
+                                                              TotalLagrangianSPHSystem},
+                                       system_coords)
+    (; smoothing_length, smoothing_length_factor) = particle_system.cache
+
+    neighborhood_search = get_neighborhood_search(particle_system, neighbor_system, semi)
+    u_neighbor_system = wrap_u(u_ode, neighbor_system, semi)
+    neighbor_coords = current_coordinates(u_neighbor_system, neighbor_system)
+
+    # Loop over all pairs of particles and neighbors within the kernel cutoff.
+    foreach_point_neighbor(particle_system, neighbor_system,
+                           system_coords, neighbor_coords,
+                           neighborhood_search) do particle, neighbor, pos_diff, distance
+        # Only consider particles with a distance > 0.
+        distance < sqrt(eps()) && return
+
+        dp_neighbor = particle_spacing(neighbor_system, neighbor)
+        dp_particle = smoothing_length[particle] / smoothing_length_factor
+
+        smoothing_length[particle] = smoothing_length_factor * min(dp_neighbor, dp_particle)
+    end
+
+    return particle_system
+end