Make kernel_grad threshold compatible with single precision and small particle spacing (#913)

LasNikas · LasNikas · web-flow · commit a5972b00b64f · 2025-09-23T10:13:54.000+02:00
* adapt threshold kernel gradient

* add link

* apply formatter

* fix tests

* add check in ddt

* consistent comments

* fix

* make volume check consistent

* fix again

* fix tests

* implement suggestions

* fix gpu tests

* fix sqrt

* fix sqrt bugs

* fix numerical precision

* implement suggestions

* fix tests

* implement suggestions

---------

Co-authored-by: LasNikas &lt;niklas.nehe@web.de&gt;
diff --git a/examples/n_body/n_body_system.jl b/examples/n_body/n_body_system.jl
@@ -59,8 +59,8 @@ function TrixiParticles.interact!(dv, v_particle_system, u_particle_system,
     TrixiParticles.foreach_point_neighbor(particle_system, neighbor_system,
                                           system_coords, neighbor_coords,
                                           semi) do particle, neighbor, pos_diff, distance
-        # Only consider particles with a distance > 0
-        distance < sqrt(eps()) && return
+        # No interaction of a particle with itself
+        particle_system === neighbor_system && particle === neighbor && return
 
         # Original version
         # dv = -G * mass[neighbor] * pos_diff / norm(pos_diff)^3
diff --git a/src/general/corrections.jl b/src/general/corrections.jl
@@ -216,7 +216,9 @@ function compute_correction_values!(system,
                        smoothing_length(system, particle))
 
             kernel_correction_coefficient[particle] += volume * W
-            if distance > sqrt(eps())
+
+            # Only consider particles with a distance > 0. See `src/general/smoothing_kernels.jl` for more details.
+            if distance^2 > eps(initial_smoothing_length(system)^2)
                 grad_W = kernel_grad(system_smoothing_kernel(system), pos_diff, distance,
                                      smoothing_length(system, particle))
                 tmp = volume * grad_W
diff --git a/src/general/smoothing_kernels.jl b/src/general/smoothing_kernels.jl
@@ -3,8 +3,13 @@ abstract type AbstractSmoothingKernel{NDIMS} end
 @inline Base.ndims(::AbstractSmoothingKernel{NDIMS}) where {NDIMS} = NDIMS
 
 @inline function kernel_grad(kernel, pos_diff, distance, h)
-    # TODO Use `eps` relative to `h` to allow scaling of simulations
-    distance < sqrt(eps(typeof(h))) && return zero(pos_diff)
+    # For `distance == 0`, the analytical gradient is zero, but the code divides by zero.
+    # To account for rounding errors, we check if `distance` is almost zero.
+    # Since the coordinates are in the order of the smoothing length `h`,
+    # `distance^2` is in the order of `h^2`, hence the comparison `distance^2 < eps(h^2)`.
+    # Note that this is faster than `distance < sqrt(eps(h^2))`.
+    # Also note that `sqrt(eps(h^2)) != eps(h)`.
+    distance^2 < eps(h^2) && return zero(pos_diff)
 
     return kernel_deriv(kernel, distance, h) / distance * pos_diff
 end
diff --git a/src/preprocessing/particle_packing/rhs.jl b/src/preprocessing/particle_packing/rhs.jl
@@ -8,8 +8,8 @@ function interact!(dv, v_particle_system, u_particle_system,
     # Loop over all pairs of particles and neighbors within the kernel cutoff
     foreach_point_neighbor(system, neighbor_system, system_coords, neighbor_coords,
                            semi) do particle, neighbor, pos_diff, distance
-        # Only consider particles with a distance > 0
-        distance < sqrt(eps()) && return
+        # Only consider particles with a distance > 0. See `src/general/smoothing_kernels.jl` for more details.
+        distance^2 < eps(initial_smoothing_length(system)^2) && return
 
         rho_a = system.initial_condition.density[particle]
         rho_b = neighbor_system.initial_condition.density[neighbor]
diff --git a/src/schemes/boundary/open_boundary/method_of_characteristics.jl b/src/schemes/boundary/open_boundary/method_of_characteristics.jl
@@ -102,7 +102,7 @@ end
 # J2: Propagates downstream to the local flow
 # J3: Propagates upstream to the local flow
 function evaluate_characteristics!(system, v, u, v_ode, u_ode, semi, t)
-    (; volume, cache, fluid_system, density, pressure) = system
+    (; volume, cache, fluid_system, smoothing_length) = system
     (; characteristics, previous_characteristics) = cache
 
     @threaded semi for particle in eachparticle(system)
@@ -165,8 +165,9 @@ function evaluate_characteristics!(system, v, u, v_ode, u_ode, semi, t)
     # of fluid particles by using the average of the values of the previous time step.
     # See eq. 27 in Negi (2020) https://doi.org/10.1016/j.cma.2020.113119
     @threaded semi for particle in each_integrated_particle(system)
-        # Particle is outside of the influence of fluid particles
-        if isapprox(volume[particle], 0)
+        # Particle is outside of the influence of fluid particles.
+        # `volume` is in the order of 1 / h^d, so volume * h^d is in the order of 1.
+        if volume[particle] * smoothing_length^ndims(system) < eps(eltype(smoothing_length))
 
             # Using the average of the values at the previous time step for particles which
             # are outside of the influence of fluid particles.
@@ -178,7 +179,9 @@ function evaluate_characteristics!(system, v, u, v_ode, u_ode, semi, t)
             for neighbor in each_integrated_particle(system)
                 # Make sure that only neighbors in the influence of
                 # the fluid particles are used.
-                if volume[neighbor] > sqrt(eps())
+                # `volume` is in the order of 1 / h^d, so volume * h^d is in the order of 1.
+                if volume[neighbor] * smoothing_length^ndims(system) >
+                   eps(eltype(smoothing_length))
                     avg_J1 += previous_characteristics[1, neighbor]
                     avg_J2 += previous_characteristics[2, neighbor]
                     avg_J3 += previous_characteristics[3, neighbor]
diff --git a/src/schemes/fluid/entropically_damped_sph/rhs.jl b/src/schemes/fluid/entropically_damped_sph/rhs.jl
@@ -67,10 +67,11 @@ function interact!(dv, v_particle_system, u_particle_system,
         dv_adhesion = adhesion_force(surface_tension_a, particle_system, neighbor_system,
                                      particle, neighbor, pos_diff, distance)
 
+        dv_particle = dv_pressure + dv_viscosity_ + dv_tvf + dv_surface_tension +
+                      dv_adhesion
+
         for i in 1:ndims(particle_system)
-            @inbounds dv[i,
-                         particle] += dv_pressure[i] + dv_viscosity_[i] + dv_tvf[i] +
-                                      dv_surface_tension[i] + dv_adhesion[i]
+            @inbounds dv[i, particle] += dv_particle[i]
         end
 
         v_diff = current_velocity(v_particle_system, particle_system, particle) -
diff --git a/src/schemes/fluid/surface_tension.jl b/src/schemes/fluid/surface_tension.jl
@@ -174,8 +174,8 @@ end
                                        neighbor_system::AbstractFluidSystem,
                                        particle, neighbor, pos_diff, distance,
                                        rho_a, rho_b, grad_kernel)
-    # No cohesion with oneself
-    distance < sqrt(eps()) && return zero(pos_diff)
+    # No cohesion with oneself. See `src/general/smoothing_kernels.jl` for more details.
+    distance^2 < eps(initial_smoothing_length(particle_system)^2) && return zero(pos_diff)
 
     m_b = hydrodynamic_mass(neighbor_system, neighbor)
     support_radius = compact_support(smoothing_kernel,
@@ -194,8 +194,8 @@ end
     (; surface_tension_coefficient) = surface_tension_a
 
     smoothing_length_ = smoothing_length(particle_system, particle)
-    # No surface tension with oneself
-    distance < sqrt(eps()) && return zero(pos_diff)
+    # No surface tension with oneself. See `src/general/smoothing_kernels.jl` for more details.
+    distance^2 < eps(initial_smoothing_length(particle_system)^2) && return zero(pos_diff)
 
     m_b = hydrodynamic_mass(neighbor_system, neighbor)
     n_a = surface_normal(particle_system, particle)
@@ -215,8 +215,8 @@ end
                                        rho_a, rho_b, grad_kernel)
     (; surface_tension_coefficient) = surface_tension_a
 
-    # No surface tension with oneself
-    distance < sqrt(eps()) && return zero(pos_diff)
+    # No surface tension with oneself. See `src/general/smoothing_kernels.jl` for more details.
+    distance^2 < eps(initial_smoothing_length(particle_system)^2) && return zero(pos_diff)
 
     n_a = surface_normal(particle_system, particle)
     curvature_a = curvature(particle_system, particle)
@@ -285,8 +285,8 @@ end
                                        rho_a, rho_b, grad_kernel)
     (; surface_tension_coefficient) = surface_tension_a
 
-    # No surface tension with oneself
-    distance < sqrt(eps()) && return zero(pos_diff)
+    # No surface tension with oneself. See `src/general/smoothing_kernels.jl` for more details.
+    distance^2 < eps(initial_smoothing_length(particle_system)^2) && return zero(pos_diff)
 
     S_a = stress_tensor(particle_system, particle)
     S_b = stress_tensor(neighbor_system, neighbor)
@@ -302,8 +302,8 @@ end
                                 pos_diff, distance)
     (; adhesion_coefficient) = neighbor_system
 
-    # No adhesion with oneself
-    distance < sqrt(eps()) && return zero(pos_diff)
+    # No adhesion with oneself. See `src/general/smoothing_kernels.jl` for more details.
+    distance^2 < eps(initial_smoothing_length(particle_system)^2) && return zero(pos_diff)
 
     # No reason to calculate the adhesion force if adhesion coefficient is near zero
     abs(adhesion_coefficient) < eps() && return zero(pos_diff)
diff --git a/src/schemes/fluid/weakly_compressible_sph/density_diffusion.jl b/src/schemes/fluid/weakly_compressible_sph/density_diffusion.jl
@@ -187,8 +187,9 @@ function update!(density_diffusion::DensityDiffusionAntuono, v, u, system, semi)
     foreach_point_neighbor(system, system, system_coords, system_coords, semi;
                            points=each_integrated_particle(system)) do particle, neighbor,
                                                                        pos_diff, distance
-        # Only consider particles with a distance > 0
-        distance < sqrt(eps(typeof(distance))) && return
+        # Density diffusion terms are all zero for distance zero.
+        # See `src/general/smoothing_kernels.jl` for more details.
+        distance^2 < eps(initial_smoothing_length(system)^2) && return
 
         rho_a = current_density(v, system, particle)
         rho_b = current_density(v, system, neighbor)
@@ -215,8 +216,9 @@ end
                                                 pos_diff, distance, m_b, rho_a, rho_b,
                                                 particle_system::AbstractFluidSystem,
                                                 grad_kernel)
-    # Density diffusion terms are all zero for distance zero
-    distance < sqrt(eps(typeof(distance))) && return
+    # Density diffusion terms are all zero for distance zero.
+    # See `src/general/smoothing_kernels.jl` for more details.
+    distance^2 < eps(initial_smoothing_length(particle_system)^2) && return
 
     (; delta) = density_diffusion
     sound_speed = system_sound_speed(particle_system)
diff --git a/src/schemes/fluid/weakly_compressible_sph/rhs.jl b/src/schemes/fluid/weakly_compressible_sph/rhs.jl
@@ -86,10 +86,11 @@ function interact!(dv, v_particle_system, u_particle_system,
         dv_adhesion = adhesion_force(surface_tension_a, particle_system, neighbor_system,
                                      particle, neighbor, pos_diff, distance)
 
+        dv_particle = dv_pressure + dv_viscosity_ + dv_tvf + dv_surface_tension +
+                      dv_adhesion
+
         for i in 1:ndims(particle_system)
-            @inbounds dv[i,
-                         particle] += dv_pressure[i] + dv_viscosity_[i] + dv_tvf[i] +
-                                      dv_surface_tension[i] + dv_adhesion[i]
+            @inbounds dv[i, particle] += dv_particle[i]
             # Debug example
             # debug_array[i, particle] += dv_pressure[i]
         end
diff --git a/src/schemes/structure/discrete_element_method/rhs.jl b/src/schemes/structure/discrete_element_method/rhs.jl
@@ -12,7 +12,8 @@ function interact!(dv, v_particle_system, u_particle_system, v_neighbor_system,
                                                                                 neighbor,
                                                                                 pos_diff,
                                                                                 distance
-        distance < sqrt(eps()) && return
+        # See `src/general/smoothing_kernels.jl` for more details.
+        distance^2 < eps(first(particle_system.radius)^2) && return
 
         # Retrieve particle properties
         m_a = particle_system.mass[particle]
diff --git a/src/schemes/structure/total_lagrangian_sph/rhs.jl b/src/schemes/structure/total_lagrangian_sph/rhs.jl
@@ -26,8 +26,8 @@ end
                            points=each_integrated_particle(system)) do particle, neighbor,
                                                                        initial_pos_diff,
                                                                        initial_distance
-        # Only consider particles with a distance > 0
-        initial_distance < sqrt(eps()) && return
+        # Only consider particles with a distance > 0. See `src/general/smoothing_kernels.jl` for more details.
+        initial_distance^2 < eps(initial_smoothing_length(system)^2) && return
 
         rho_a = @inbounds system.material_density[particle]
         rho_b = @inbounds system.material_density[neighbor]
@@ -55,10 +55,10 @@ end
                                           current_pos_diff, current_distance,
                                           m_a, m_b, rho_a, rho_b, grad_kernel)
 
+        dv_particle = dv_stress + dv_penalty_force_ + dv_viscosity
+
         for i in 1:ndims(system)
-            @inbounds dv[i,
-                         particle] += dv_stress[i] + dv_penalty_force_[i] +
-                                      dv_viscosity[i]
+            @inbounds dv[i, particle] += dv_particle[i]
         end
 
         # TODO continuity equation for boundary model with `ContinuityDensity`?
@@ -88,8 +88,8 @@ function interact!(dv, v_particle_system, u_particle_system,
                                                                                 neighbor,
                                                                                 pos_diff,
                                                                                 distance
-        # Only consider particles with a distance > 0
-        distance < sqrt(eps()) && return
+        # Only consider particles with a distance > 0. See `src/general/smoothing_kernels.jl` for more details.
+        distance^2 < eps(initial_smoothing_length(particle_system)^2) && return
 
         # Apply the same force to the structure particle
         # that the fluid particle experiences due to the structure particle.
diff --git a/src/schemes/structure/total_lagrangian_sph/system.jl b/src/schemes/structure/total_lagrangian_sph/system.jl
@@ -326,8 +326,8 @@ end
     initial_coords = initial_coordinates(system)
     foreach_point_neighbor(system, system, initial_coords, initial_coords,
                            semi) do particle, neighbor, initial_pos_diff, initial_distance
-        # Only consider particles with a distance > 0.
-        initial_distance < sqrt(eps()) && return
+        # Only consider particles with a distance > 0. See `src/general/smoothing_kernels.jl` for more details.
+        initial_distance^2 < eps(initial_smoothing_length(system)^2) && return
 
         volume = mass[neighbor] / material_density[neighbor]
         pos_diff = current_coords(system, particle) - current_coords(system, neighbor)
diff --git a/test/schemes/structure/total_lagrangian_sph/rhs.jl b/test/schemes/structure/total_lagrangian_sph/rhs.jl
@@ -85,13 +85,13 @@
 
             TrixiParticles.eachparticle(::MockSystem) = eachparticle
             TrixiParticles.each_integrated_particle(::MockSystem) = each_integrated_particle
+            TrixiParticles.smoothing_length(::MockSystem, _) = eps()
 
             function TrixiParticles.add_acceleration!(_, _, ::MockSystem)
                 return nothing
             end
             TrixiParticles.kernel_deriv(::Val{:mock_smoothing_kernel}, _, _) = kernel_deriv
             TrixiParticles.compact_support(::MockSystem, ::MockSystem) = 1000.0
-            Base.eps(::Type{Val{:mock_smoothing_length}}) = eps()
             function TrixiParticles.current_coords(system::MockSystem, particle)
                 return TrixiParticles.current_coords(initial_coordinates, system, particle)
             end
