Clarify and visualize SPH surface normals

NEWS.md

@@ -4,6 +4,14 @@ TrixiParticles.jl follows the interpretation of
 [semantic versioning (semver)](https://julialang.github.io/Pkg.jl/dev/compatibility/#Version-specifier-format-1)
 used in the Julia ecosystem. Notable changes will be documented in this file for human readability.
 
+## Version 0.5.4
+
+### Important Bugfixes
+
+- Hardened surface tension model configuration by validating coefficients and surface-normal
+  thresholds, avoiding unnecessary normal allocation for `CohesionForceAkinci`, and stabilizing
+  Akinci cohesion and adhesion kernels across floating-point scales.
+
 ## Version 0.5.3
 
 ### Features

docs/src/systems/fluid.md

@@ -215,17 +215,20 @@ Pages = [joinpath("general", "corrections.jl")]
 
 ### Overview of surface normal calculation in SPH
 
-Surface normals are essential for modeling surface tension as they provide the directionality
-of forces acting at the fluid interface. They are calculated based on the particle properties and
-their spatial distribution.
+Surface normals characterize the local orientation of an interface. In SPH, this geometric
+information can be used for interface detection and reconstruction, curvature estimation,
+interfacial boundary conditions, and interfacial force models. The computed normal field is also
+available for analysis and VTK output.
 
 #### Color field and gradient-based surface normals
 
-The surface normal at a particle is derived from the color field, a scalar field assigned to particles
-to distinguish between different fluid phases or between fluid and air. The color field gradients point
-towards the interface, and the normalized gradient defines the surface normal direction.
+The surface normal at a particle can be derived from a color field, a scalar marker used to
+distinguish phases or materials. Its gradient is perpendicular to the color-field level sets and
+therefore provides an interface-normal estimate; its orientation depends on the chosen color
+convention. For a free surface whose exterior phase is not represented by particles, truncation of
+the kernel support creates the corresponding discrete color-field gradient.
 
-The simplest SPH formulation for a surface normal, ``n_a`` is given as
+The simplest SPH approximation of an unnormalized color-field normal, ``n_a``, is
 
 ```math
 n_a = \sum_b m_b \frac{c_b}{\rho_b} \nabla_a W_{ab},
@@ -238,23 +241,117 @@ where:
 - ``\rho_b`` is the density of particle ``b``,
 - ``\nabla_a W_{ab}`` is the gradient of the smoothing kernel ``W_{ab}`` with respect to particle ``a``.
 
+```@eval
+using CairoMakie
+
+let
+    coordinate = range(-2.0, 2.0, length=401)
+    interface_width = 0.3
+    colorfield = @. 0.5 * (1.0 - tanh(coordinate / interface_width))
+    colorfield_gradient = @. -0.5 / interface_width /
+                             cosh(coordinate / interface_width)^2
+
+    fig = Figure(size=(1000, 430), fontsize=18)
+    color_axis = Axis(fig[1, 1],
+                      xlabel="signed distance s/h", ylabel="color field c",
+                      title="Diffuse color-field transition")
+    gradient_axis = Axis(fig[1, 2],
+                         xlabel="signed distance s/h", ylabel="dc/d(s/h)",
+                         title="Color-field gradient")
+
+    lines!(color_axis, coordinate, colorfield, color=:steelblue, linewidth=3)
+    lines!(gradient_axis, coordinate, colorfield_gradient, color=:darkorange,
+           linewidth=3)
+    vlines!(color_axis, [0.0], color=:black, linestyle=:dash, linewidth=2)
+    vlines!(gradient_axis, [0.0], color=:black, linestyle=:dash, linewidth=2)
+    hlines!(gradient_axis, [0.0], color=(:black, 0.35), linewidth=1)
+    xlims!(color_axis, extrema(coordinate))
+    xlims!(gradient_axis, extrema(coordinate))
+
+    CairoMakie.save("colorfield_profile.png", fig)
+end
+```
+
+![A diffuse color field and its gradient across an interface](colorfield_profile.png)
+
+The color field is approximately constant within either phase. Its gradient is localized in the
+transition region and vanishes away from the interface. The sign of the gradient determines the
+normal orientation; exchanging the two color values reverses that orientation.
+
 #### Normalization of surface normals
 
-The calculated normals are normalized to unit vectors:
+The color-field gradient ``n_a`` is generally not a unit vector. Formulations that require only
+the interface orientation use the unit normal
 
 ```math
 \hat{n}_a = \frac{n_a}{\Vert n_a \Vert}.
 ```
 
-Normalization ensures that the magnitude of the normals does not bias the curvature calculations or the resulting surface tension forces.
+Normalization separates the interface orientation from the magnitude of the discrete color-field
+gradient. In TrixiParticles, the Morris formulations use unit normals for curvature or surface-stress
+calculations, while the Akinci surface-area force uses the unnormalized gradient.
+
+```@eval
+using CairoMakie
+
+let
+    coordinate = range(-1.35, 1.35, length=241)
+    radius = 0.85
+    interface_width = 0.1
+    colorfield = [0.5 * (1.0 - tanh((hypot(x, y) - radius) / interface_width))
+                  for x in coordinate, y in coordinate]
+
+    particle_spacing = 0.17
+    particle_coordinates = [(x, y) for x in (-radius):particle_spacing:radius
+                             for y in (-radius):particle_spacing:radius
+                             if hypot(x, y) <= radius]
+    particle_x = first.(particle_coordinates)
+    particle_y = last.(particle_coordinates)
+
+    angles = range(0.0, 2pi, length=13)[1:(end - 1)]
+    normal_x = -cos.(angles)
+    normal_y = -sin.(angles)
+    interface_x = radius .* cos.(angles)
+    interface_y = radius .* sin.(angles)
+
+    fig = Figure(size=(760, 650), fontsize=18)
+    axis = Axis(fig[1, 1], aspect=DataAspect(),
+                xlabel="x/h", ylabel="y/h",
+                title="Interface orientation from the color-field gradient")
+    heatmap = heatmap!(axis, coordinate, coordinate, colorfield,
+                       colormap=:viridis, colorrange=(0.0, 1.0))
+    contour!(axis, coordinate, coordinate, colorfield,
+             levels=[0.1, 0.9], color=(:white, 0.8), linewidth=1.5)
+    contour!(axis, coordinate, coordinate, colorfield,
+             levels=[0.5], color=:white, linewidth=3)
+    scatter!(axis, particle_x, particle_y, color=(:black, 0.35), markersize=5)
+    arrows2d!(axis, interface_x, interface_y, normal_x, normal_y,
+              normalize=true, lengthscale=0.3, color=:black,
+              shaftwidth=3, tipwidth=14, tiplength=10)
+    xlims!(axis, extrema(coordinate))
+    ylims!(axis, extrema(coordinate))
+    Colorbar(fig[1, 2], heatmap, label="color field c")
+
+    CairoMakie.save("colorfield_surface_normals.png", fig)
+end
+```
+
+![Color-field level sets and interface-normal directions](colorfield_surface_normals.png)
+
+The particle phase has ``c \approx 1`` and the exterior has ``c \approx 0``. Consequently,
+``\nabla c`` and the displayed unit normals point toward increasing ``c``. The arrows are
+perpendicular to the color-field level sets; reversing the color convention reverses the arrows
+without changing the interface geometry.
 
 #### Handling noise and errors in normal calculation
 
-In regions distant from the interface, the calculated normals may be small or inaccurate due to the
-smoothing kernel's support radius. To mitigate this:
+Away from an interface, the exact color-field gradient vanishes, but particle disorder and
+incomplete kernel support can produce small or poorly resolved normal estimates. To mitigate this:
 
-1. Normals below a threshold are excluded from further calculations.
-2. Curvature calculations use a corrected formulation to reduce errors near interface fringes.
+1. Normals with insufficient particle support are discarded.
+2. Morris-type calculations reject small gradients and can suppress normals in well-resolved
+   interior regions.
+3. Curvature calculations use a corrected formulation to reduce errors near interface fringes.
 
 ```@autodocs
 Modules = [TrixiParticles]
@@ -289,6 +386,17 @@ In the following table some values are shown for reference. The values marked wi
 | **Water**       | 0.07288  [Lange](@cite Lange2005)               |
 | **Mercury**     | 0.486502 [Lange](@cite Lange2005)               |
 
+### Model configuration
+
+All surface tension coefficients must be finite and non-negative. A zero coefficient disables
+the fluid-fluid surface force. Wall adhesion is controlled independently by the boundary's
+`adhesion_coefficient`.
+
+`CohesionForceAkinci` only evaluates the pairwise cohesion and optional wall-adhesion forces.
+It does not require surface normals or `reference_particle_spacing`. The full
+`SurfaceTensionAkinci` model and both Morris models require a surface-normal method. When one
+of these models is selected without an explicit method, `ColorfieldSurfaceNormal()` is used.
+
 ### [Akinci-based intra-particle force surface tension and wall adhesion model](@id akinci_ipf)
 
 The [Akinci](@cite Akinci2013) model divides surface tension into distinct force components:
@@ -350,6 +458,10 @@ A(r) = \frac{0.007}{h_c^{3.25}}
 \end{cases}
 ```
 
+The published adhesion kernel uses a three-dimensional normalization. In two-dimensional
+simulations, `adhesion_coefficient` is therefore an empirical numerical parameter and may need
+to be adjusted when changing the particle spacing or smoothing length.
+
 ---
 
 ### [Morris surface tension model](@id morris_csf)

src/general/semidiscretization.jl

@@ -920,8 +920,7 @@ function check_system_color(systems)
         system isa AbstractFluidSystem || return false
         system isa ParticlePackingSystem && return false
 
-        return !isnothing(system.surface_tension) ||
-               system.surface_normal_method isa ColorfieldSurfaceNormal
+        return system.surface_normal_method isa ColorfieldSurfaceNormal
     end
 
     if requires_color_check
@@ -935,7 +934,7 @@ function check_system_color(systems)
                              systems)
 
         if length(system_ids) > 1 && sum(i -> systems[i].cache.color, system_ids) == 0
-            throw(ArgumentError("If `ColorfieldSurfaceNormal` or a surface tension model is used, at least one participating system must have a color different from 0."))
+            throw(ArgumentError("When `ColorfieldSurfaceNormal` is used, at least one participating system must have a color different from 0."))
         end
     end
 end

src/schemes/boundary/wall_boundary/dummy_particles.jl

@@ -33,7 +33,7 @@ Boundary model for [`WallBoundarySystem`](@ref).
                                 in areas of low pressure, against which the particle
                                 shifting technique is fighting.
 - `reference_particle_spacing`: The reference particle spacing used for weighting values at the boundary,
-                                which currently is only needed when using surface tension.
+                                which is needed when using a surface-normal method.
 # Examples
 ```jldoctest; output = false, setup = :(densities = [1.0, 2.0, 3.0]; masses = [0.1, 0.2, 0.3]; smoothing_kernel = SchoenbergCubicSplineKernel{2}(); smoothing_length = 0.1)
 # Free-slip condition

src/schemes/fluid/entropically_damped_sph/system.jl

@@ -50,10 +50,13 @@ See [Entropically Damped Artificial Compressibility for SPH](@ref edac) for more
                                 The keyword argument `acceleration` should be used instead for
                                 gravity-like source terms.
 - `surface_tension`:            Surface tension model used for this SPH system. (default: no surface tension)
-- `surface_normal_method`:      The surface normal method to be used for this SPH system.
-                                (default: no surface normal method or `ColorfieldSurfaceNormal()` if a surface_tension model is used)
+- `surface_normal_method`:      Method used to estimate fluid-interface normals. This can be
+                                configured independently for interface analysis and output and is
+                                also used by models that require interface geometry. The default is
+                                `nothing`; `ColorfieldSurfaceNormal()` is selected automatically
+                                when the surface tension model requires normals.
 - `reference_particle_spacing`: The reference particle spacing used for weighting values at the boundary,
-                                which currently is only needed when using surface tension.
+                                which is needed when using a surface-normal method.
 - `color_value`:                Integer label used for calculation of surface normals.
                                 Currently this is only used together with [`BoundaryModelDummyParticles`](@ref) and
                                 [`ColorfieldSurfaceNormal`](@ref): fluid-boundary normal evaluation
@@ -119,12 +122,11 @@ function EntropicallyDampedSPHSystem(initial_condition; smoothing_kernel, smooth
         throw(ArgumentError("`acceleration` must be of length $NDIMS for a $(NDIMS)D problem"))
     end
 
-    if surface_tension !== nothing && surface_normal_method === nothing
-        surface_normal_method = ColorfieldSurfaceNormal()
-    end
+    surface_normal_method = default_surface_normal_method(surface_tension,
+                                                          surface_normal_method)
 
     if surface_normal_method !== nothing && reference_particle_spacing < eps()
-        throw(ArgumentError("`reference_particle_spacing` must be set to a positive value when using `ColorfieldSurfaceNormal` or a surface tension model"))
+        throw(ArgumentError("`reference_particle_spacing` must be set to a positive value when using a surface-normal method"))
     end
 
     if correction isa ShepardKernelCorrection &&
@@ -308,7 +310,7 @@ function update_final!(system::EntropicallyDampedSPHSystem, v, u, v_ode, u_ode,
                        kwargs...)
     (; surface_tension) = system
 
-    # Surface normal of neighbor and boundary needs to have been calculated already
+    # Surface-tension formulations using curvature or stress require previously computed normals.
     compute_curvature!(system, surface_tension, v, u, v_ode, u_ode, semi, t)
     compute_stress_tensors!(system, surface_tension, v, u, v_ode, u_ode, semi, t)
     update_average_pressure!(system, system.average_pressure_reduction, v_ode, u_ode, semi)

src/schemes/fluid/fluid.jl

@@ -229,10 +229,13 @@ function calculate_dt(v_ode, u_ode, cfl_number, system::AbstractFluidSystem, sem
 
     if surface_tension isa SurfaceTensionMorris ||
        surface_tension isa SurfaceTensionMomentumMorris
-        v = wrap_v(v_ode, system, semi)
-        dt_surface_tension = sqrt(current_density(v, system, 1) * smoothing_length_^3 /
-                                  (2 * pi * surface_tension.surface_tension_coefficient))
-        dt = min(dt, dt_surface_tension)
+        coefficient = surface_tension.surface_tension_coefficient
+        if !iszero(coefficient)
+            v = wrap_v(v_ode, system, semi)
+            dt_surface_tension = sqrt(current_density(v, system, 1) * smoothing_length_^3 /
+                                      (2 * pi * coefficient))
+            dt = min(dt, dt_surface_tension)
+        end
     end
 
     return dt

src/schemes/fluid/surface_normal_sph.jl

@@ -2,7 +2,9 @@
     ColorfieldSurfaceNormal(; boundary_contact_threshold=0.1, interface_threshold=0.01,
                               ideal_density_threshold=0.0)
 
-Color field based computation of the interface normals.
+Color-field-based computation of fluid-interface normals. Interface normals describe local
+interface geometry and can be computed for analysis, output, or use by models that require
+interface orientation.
 
 # Keywords
 - `boundary_contact_threshold=0.1`: If this threshold is reached the fluid is assumed to be in contact with the boundary.
@@ -17,8 +19,22 @@ end
 
 function ColorfieldSurfaceNormal(; boundary_contact_threshold=0.1, interface_threshold=0.01,
                                  ideal_density_threshold=0.0)
-    return ColorfieldSurfaceNormal(boundary_contact_threshold, interface_threshold,
-                                   ideal_density_threshold)
+    thresholds = (boundary_contact_threshold, interface_threshold,
+                  ideal_density_threshold)
+    if !all(threshold -> threshold isa Real && isfinite(threshold), thresholds)
+        throw(ArgumentError("surface-normal thresholds must be finite real numbers"))
+    end
+
+    thresholds = promote(thresholds...)
+    return ColorfieldSurfaceNormal(thresholds...)
+end
+
+@inline function default_surface_normal_method(surface_tension, surface_normal_method)
+    if isnothing(surface_normal_method) && requires_surface_normal(surface_tension)
+        return ColorfieldSurfaceNormal()
+    end
+
+    return surface_normal_method
 end
 
 function create_cache_surface_normal(surface_normal_method, ELTYPE, NDIMS, nparticles)