fix

src/TrixiParticles.jl

@@ -41,6 +41,7 @@ include("general/system.jl")
 # `util.jl` needs to be next because of the macros `@trixi_timeit` and `@threaded`
 include("util.jl")
 include("preprocessing/preprocessing.jl")
+include("multi_resolution/multi_resolution.jl")
 include("callbacks/callbacks.jl")
 include("general/general.jl")
 include("setups/setups.jl")

src/multi_resolution/multi_resolution.jl

@@ -1,3 +1,4 @@
+include("refinement_criteria.jl")
 include("particle_refinement.jl")
 include("refinement_pattern.jl")
 

src/multi_resolution/particle_refinement.jl

@@ -19,7 +19,10 @@ function ParticleRefinement(; splitting_pattern, max_spacing_ratio,
     split_candidates = Vector{ELTYPE}()
     delete_candidates = Vector{Bool}()
 
-    return ParticleRefinement(splitting_pattern, Tuple(refinement_criteria),
+    if refinement_criteria isa Tuple
+        refinement_criteria = (refinement_criteria,)
+    end
+    return ParticleRefinement(splitting_pattern, refinement_criteria,
                               max_spacing_ratio, mass_ref, merge_candidates,
                               split_candidates, delete_candidates, 0, 0)
 end
@@ -45,12 +48,18 @@ function refinement!(semi, v_ode, u_ode, v_tmp, u_tmp, t)
 
     # Resize neighborhood search
     foreach_system(semi) do system
-        foreach_system(semi) do neighbor_system
-            search = get_neighborhood_search(system, neighbor_system, semi)
-            u_neighbor = wrap_u(u_ode, neighbor_system, semi)
+        foreach_system(semi) do neighbor
+            nhs_old = get_neighborhood_search(system, neighbor, semi)
+            u_neighbor = wrap_u(u_ode, neighbor, semi)
 
             # TODO
-            resize_nhs!(search, system, neighbor_system, u_neighbor)
+            resize_nhs!(nhs_old, system, neighbor, u_neighbor)
+            nhs = copy_neighborhood_search(nhs_old, compact_support(system, neighbor),
+                                           nparticles(neighbor))
+
+            system_coords = current_coordinates(u_ode, system)
+            neighbor_coords = current_coordinates(u_ode, neighbor)
+            PointNeighbors.initialize!(nhs, system_coords, neighbor_coords)
         end
     end
 

src/multi_resolution/refinement_pattern.jl

@@ -17,7 +17,7 @@ struct CubicSplitting{ELTYPE}
             SVector{2, ELTYPE}(epsilon * direction_3),
             SVector{2, ELTYPE}(epsilon * direction_4)]
 
-        if refinement_pattern.center_particle
+        if center_particle
             push!(relative_position, SVector(zero(ELTYPE), zero(ELTYPE)))
         end
 
@@ -45,7 +45,7 @@ struct TriangularSplitting{ELTYPE}
             SVector{2, ELTYPE}(epsilon * direction_2),
             SVector{2, ELTYPE}(epsilon * direction_3)]
 
-        if refinement_pattern.center_particle
+        if center_particle
             push!(relative_position, SVector(zero(ELTYPE), zero(ELTYPE)))
         end
 
@@ -80,7 +80,7 @@ struct HexagonalSplitting{ELTYPE}
             SVector{2, ELTYPE}(epsilon * direction_5),
             SVector{2, ELTYPE}(epsilon * direction_6)]
 
-        if refinement_pattern.center_particle
+        if center_particle
             push!(relative_position, SVector(zero(ELTYPE), zero(ELTYPE)))
         end
 

src/schemes/fluid/entropically_damped_sph/rhs.jl

@@ -157,6 +157,8 @@ end
 function pressure_damping_term(particle_system, neighbor_system, particle_refinement,
                                particle, neighbor, pos_diff, distance, beta_inv_a,
                                m_a, m_b, p_a, p_b, rho_a, rho_b)
+    (; sound_speed) = particle_system
+
     # EDAC pressure evolution
     pressure_diff = p_a - p_b
 
@@ -227,14 +229,23 @@ end
     return zero(pos_diff)
 end
 
+#TODO
 @inline function velocity_correction(particle_system, neighbor_system,
-                                     particle_refinement,
-                                     pos_diff, distance,
+                                     particle_refinement, pos_diff, distance,
+                                     v_particle_system, v_neighbor_system,
+                                     rho_a, rho_b, m_a, m_b, particle, neighbor,
+                                     grad_kernel)
+
+    return zero(grad_kernel)
+end
+
+@inline function velocity_correction(particle_system, neighbor_system::FluidSystem,
+                                     particle_refinement, pos_diff, distance,
                                      v_particle_system, v_neighbor_system,
                                      rho_a, rho_b, m_a, m_b, particle, neighbor,
                                      grad_kernel)
-    momentum_velocity_a = current_velocity(v_particle_system, system, particle)
-    advection_velocity_a = advection_velocity(v_particle_system, system, particle)
+    momentum_velocity_a = current_velocity(v_particle_system, particle_system, particle)
+    advection_velocity_a = advection_velocity(v_particle_system, particle_system, particle)
 
     momentum_velocity_b = current_velocity(v_neighbor_system, neighbor_system, neighbor)
     advection_velocity_b = advection_velocity(v_neighbor_system, neighbor_system, neighbor)

src/schemes/fluid/entropically_damped_sph/system.jl

@@ -209,16 +209,22 @@ function update_quantities!(system::EntropicallyDampedSPHSystem, v, u,
                             v_ode, u_ode, semi, t)
     compute_density!(system, u, u_ode, semi, system.density_calculator)
     update_average_pressure!(system, system.transport_velocity, v_ode, u_ode, semi)
-    factor_for_variable_smoothing_length!(system, system.particle_refinement, u_ode, semi)
+    factor_for_variable_smoothing_length!(system, system.particle_refinement,
+                                          v_ode, u_ode, semi)
 end
 
-factor_for_variable_smoothing_length!(system, ::Nothing, u_ode, semi) = system
+factor_for_variable_smoothing_length!(system, ::Nothing, v_ode, u_ode, semi) = system
 
-function factor_for_variable_smoothing_length!(system, particle_refinement, u_ode, semi)
+function factor_for_variable_smoothing_length!(system, particle_refinement,
+                                               v_ode, u_ode, semi)
     (; smoothing_kernel) = system
-    u = wrap_u(u_ode, system, semi)
+    (; beta) = system.cache
+
     NDIMS = ndims(system)
 
+    u = wrap_u(u_ode, system, semi)
+    v = wrap_v(v_ode, system, semi)
+
     # Use all other systems for the average pressure
     foreach_system(semi) do neighbor_system
         u_neighbor_system = wrap_u(u_ode, neighbor_system, semi)
@@ -238,7 +244,7 @@ function factor_for_variable_smoothing_length!(system, particle_refinement, u_od
             h_a = smoothing_length(system, particle)
             w_deriv = kernel_deriv(smoothing_kernel, distance, h_a)
 
-            beta[partcle] += -m_a * distance * w_deriv / (NDIMS * rho_a)
+            beta[particle] += -m_a * distance * w_deriv / (NDIMS * rho_a)
         end
     end
 

src/schemes/fluid/fluid.jl

@@ -19,8 +19,10 @@ end
 
 function create_cache_refinement(initial_condition, refinement, smoothing_length)
     smoothng_length_factor = smoothing_length / initial_condition.particle_spacing
+
     return (; smoothing_length=smoothing_length * ones(length(initial_condition.density)),
-            smoothing_lengt_factor=smoothng_length_factor)
+            smoothing_lengt_factor=smoothng_length_factor,
+            beta=zeros(nparticles(initial_condition)))
 end
 
 @inline hydrodynamic_mass(system::FluidSystem, particle) = system.mass[particle]

src/schemes/fluid/transport_velocity.jl

@@ -113,7 +113,7 @@ end
     momentum_velocity_b = current_velocity(v_neighbor_system, neighbor_system, neighbor)
     advection_velocity_b = advection_velocity(v_neighbor_system, neighbor_system, neighbor)
 
-    factor_a = beta_correction(particle_system, particle) / rho_a
+    factor_a = beta_correction(system, particle) / rho_a
     factor_b = beta_correction(neighbor_system, neighbor) / rho_b
 
     A_a = factor_a * momentum_velocity_a * (advection_velocity_a - momentum_velocity_a)'

src/schemes/fluid/viscosity.jl

@@ -228,7 +228,7 @@ end
     beta_inv_a = beta_correction(particle_system, particle)
 
     nu_a = kinematic_viscosity(particle_system,
-                               viscosity_model(neighbor_system, particle_system))
+                               viscosity_model(neighbor_system, particle_system), particle)
 
     grad_kernel_a = smoothing_kernel_grad(particle_system, pos_diff, distance, particle)
     grad_kernel_b = smoothing_kernel_grad(neighbor_system, pos_diff, distance, neighbor)
@@ -237,6 +237,10 @@ end
 
     tmp = (distance^2 + 0.001 * smoothing_length(particle_system, particle)^2)
 
+    v_a = viscous_velocity(v_particle_system, particle_system, particle)
+    v_b = viscous_velocity(v_neighbor_system, neighbor_system, neighbor)
+    v_diff = v_a - v_b
+
     return m_b * beta_inv_a * 4 * nu_a * dot(pos_diff, grad_W_avg) * v_diff /
            (tmp * (rho_a + rho_b))
 end
@@ -253,11 +257,9 @@ end
 
     # TODO is that correct with the smoothing lengths?
     nu_a = kinematic_viscosity(particle_system,
-                               viscosity_model(neighbor_system, particle_system),
-                               smoothing_length_particle)
+                               viscosity_model(neighbor_system, particle_system), particle)
     nu_b = kinematic_viscosity(neighbor_system,
-                               viscosity_model(particle_system, neighbor_system),
-                               smoothing_length_neighbor)
+                               viscosity_model(particle_system, neighbor_system), particle)
 
     v_a = viscous_velocity(v_particle_system, particle_system, particle)
     v_b = viscous_velocity(v_neighbor_system, neighbor_system, neighbor)