fix beta correction

src/multi_resolution/resize.jl

@@ -143,7 +143,7 @@ function deleteat!(system::FluidSystem, refinement, v, u)
             system.mass[particle] = mass_keep
             system.cache.smoothing_length[particle] = smoothing_length_keep
 
-            set_particle_pressure!( v, system,particle, pressure_keep)
+            set_particle_pressure!(v, system, particle, pressure_keep)
 
             set_particle_density!(v, system, particle, density_keep)
 

src/schemes/fluid/entropically_damped_sph/rhs.jl

@@ -64,7 +64,8 @@ function interact!(dv, v_particle_system, u_particle_system,
         v_diff = current_velocity(v_particle_system, particle_system, particle) -
                  current_velocity(v_neighbor_system, neighbor_system, neighbor)
 
-        pressure_evolution!(dv, particle_system, neighbor_system, v_diff, grad_kernel,
+        pressure_evolution!(dv, particle_system, neighbor_system,
+                            v_particle_system, v_neighbor_system, v_diff, grad_kernel,
                             particle, neighbor, pos_diff, distance, sound_speed,
                             m_a, m_b, p_a, p_b, rho_a, rho_b)
 
@@ -79,6 +80,7 @@ function interact!(dv, v_particle_system, u_particle_system,
 end
 
 @inline function pressure_evolution!(dv, particle_system, neighbor_system,
+                                     v_particle_system, v_neighbor_system,
                                      v_diff, grad_kernel, particle, neighbor,
                                      pos_diff, distance, sound_speed,
                                      m_a, m_b, p_a, p_b, rho_a, rho_b)
@@ -87,10 +89,8 @@ end
     # This is basically the continuity equation times `sound_speed^2`
     artificial_eos = m_b * rho_a / rho_b * sound_speed^2 * dot(v_diff, grad_kernel)
 
-    beta_inv_a = beta_correction(particle_system, particle_system.particle_refinement,
-                                 particle)
-    beta_inv_b = beta_correction(neighbor_system, neighbor_system.particle_refinement,
-                                 neighbor)
+    beta_inv_a = beta_correction(particle_system, particle)
+    beta_inv_b = beta_correction(neighbor_system, neighbor)
 
     dv[end, particle] += beta_inv_a * artificial_eos +
                          pressure_damping_term(particle_system, neighbor_system,
@@ -106,12 +106,18 @@ end
     return dv
 end
 
-function beta_correction(particle_system, ::Nothing, particle)
-    return one(eltype(particle_system))
+beta_correction(system, particle) = one(eltype(system))
+
+function beta_correction(system::FluidSystem, particle)
+    beta_correction(system, system.particle_refinement, particle)
 end
 
-function beta_correction(particle_system, particle_refinement, particle)
-    return inv(particle_system.cache.beta[particle])
+function beta_correction(system::FluidSystem, ::Nothing, particle)
+    return one(eltype(system))
+end
+
+function beta_correction(system, particle_refinement, particle)
+    return inv(system.cache.beta[particle])
 end
 
 function pressure_damping_term(particle_system, neighbor_system, ::Nothing,
@@ -132,6 +138,8 @@ function pressure_damping_term(particle_system, neighbor_system, ::Nothing,
                                 smoothing_length(particle_system, particle))
     tmp = eta_tilde / (distance^2 + 0.01 * smoothing_length_average^2)
 
+    grad_kernel = smoothing_kernel_grad(particle_system, pos_diff, distance, particle)
+
     # This formulation was introduced by Hu and Adams (2006). https://doi.org/10.1016/j.jcp.2005.09.001
     # They argued that the formulation is more flexible because of the possibility to formulate
     # different inter-particle averages or to assume different inter-particle distributions.
@@ -182,8 +190,8 @@ function pressure_reduction(particle_system, neighbor_system, particle_refinemen
                             v_particle_system, v_neighbor_system,
                             particle, neighbor, pos_diff, distance, m_b,
                             p_a, p_b, rho_a, rho_b, beta_a, beta_b)
-    beta_inv_a = beta_correction(particle_system, particle_refinement, particle)
-    beta_inv_b = beta_correction(particle_system, particle_refinement, neighbor)
+    beta_inv_a = beta_correction(particle_system, particle)
+    beta_inv_b = beta_correction(particle_system, neighbor)
 
     p_a_avg = average_pressure(particle_system, particle)
     p_b_avg = average_pressure(neighbor_system, neighbor)
@@ -234,7 +242,7 @@ end
     v_diff = momentum_velocity_a - momentum_velocity_b
     v_diff_tilde = advection_velocity_a - advection_velocity_b
 
-    beta_inv_a = beta_correction(particle_system, particle_refinement, particle)
+    beta_inv_a = beta_correction(particle_system, particle)
 
     return -m_b * beta_inv_a *
            (dot(v_diff_tilde - v_diff, grad_kernel) * momentum_velocity_a) / rho_b
@@ -285,10 +293,8 @@ end
     p_a_avg = average_pressure(particle_system, particle)
     p_b_avg = average_pressure(neighbor_system, neighbor)
 
-    P_a = beta_correction(particle_system, particle_refinement, particle) *
-          (p_a - p_a_avg) / rho_a^2
-    P_b = beta_correction(particle_system, particle_refinement, neighbor) *
-          (p_b - p_b_avg) / rho_b^2
+    P_a = beta_correction(particle_system, particle) * (p_a - p_a_avg) / rho_a^2
+    P_b = beta_correction(particle_system, neighbor) * (p_b - p_b_avg) / rho_b^2
 
     grad_kernel_a = smoothing_kernel_grad(particle_system, pos_diff, distance, particle)
     grad_kernel_b = smoothing_kernel_grad(neighbor_system, pos_diff, distance, neighbor)

src/schemes/fluid/transport_velocity.jl

@@ -113,8 +113,8 @@ 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_refinement, particle) / rho_a
-    factor_b = beta_correction(neighbor_system, particle_refinement, neighbor) / rho_b
+    factor_a = beta_correction(particle_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)'
     A_b = factor_b * momentum_velocity_b * (advection_velocity_b - momentum_velocity_b)'

src/schemes/fluid/viscosity.jl

@@ -225,7 +225,7 @@ end
                                              particle, neighbor, pos_diff,
                                              distance, sound_speed, m_a, m_b,
                                              rho_a, rho_b, grad_kernel)
-    beta_inv_a = beta_correction(particle_system, particle_refinement, particle)
+    beta_inv_a = beta_correction(particle_system, particle)
 
     nu_a = kinematic_viscosity(particle_system,
                                viscosity_model(neighbor_system, particle_system))