Move mean density calculation into free surface calculation

Author: efaulhaber

src/general/corrections.jl

@@ -35,15 +35,16 @@ end
 # `rho_mean` is the mean density of the fluid, which is used to determine correction values near the free surface.
 #  Return a tuple `(viscosity_correction, pressure_correction, surface_tension_correction)` representing the correction terms.
 @inline function free_surface_correction(correction::AkinciFreeSurfaceCorrection,
-                                         particle_system, rho_mean)
+                                         particle_system, rho_a, rho_b)
     # Equation 4 in ref
+    rho_mean = (rho_a + rho_b) / 2
     k = correction.rho0 / rho_mean
 
     # Viscosity, pressure, surface_tension
     return k, 1, k
 end
 
-@inline function free_surface_correction(correction, particle_system, rho_mean)
+@inline function free_surface_correction(correction, particle_system, rho_a, rho_b)
     return 1, 1, 1
 end
 

src/schemes/fluid/weakly_compressible_sph/rhs.jl

@@ -58,7 +58,6 @@ function interact!(dv, v_particle_system, u_particle_system,
             m_b = @inbounds hydrodynamic_mass(neighbor_system, neighbor)
             v_b = @inbounds current_velocity(v_neighbor_system, neighbor_system, neighbor)
             rho_b = @inbounds current_density(v_neighbor_system, neighbor_system, neighbor)
-            rho_mean = (rho_a + rho_b) / 2
 
             # The following call is equivalent to
             #     `p_b = current_pressure(v_neighbor_system, neighbor_system, neighbor)`
@@ -73,7 +72,7 @@ function interact!(dv, v_particle_system, u_particle_system,
             (viscosity_correction, pressure_correction,
              surface_tension_correction) = free_surface_correction(correction,
                                                                    particle_system,
-                                                                   rho_mean)
+                                                                   rho_a, rho_b)
 
             # For `ContinuityDensity` without correction, this is equivalent to
             # dv_pressure = -m_b * (p_a + p_b) / (rho_a * rho_b) * grad_kernel