Implement multicolor surface normals and threshold handling

NEWS.md

@@ -4,6 +4,16 @@ 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.
+- Added color-weighted fluid-fluid interface normals and consistently applied surface-normal
+  validity thresholds to standalone, Akinci, and Morris calculations.
+
 ## Version 0.5.3
 
 ### Features

src/general/semidiscretization.jl

@@ -917,25 +917,29 @@ check_configuration(system::AbstractSystem, systems, nhs) = nothing
 
 function check_system_color(systems)
     requires_color_check = any(systems) do system
-        system isa AbstractFluidSystem || return false
-        system isa ParticlePackingSystem && return false
-
-        return !isnothing(system.surface_tension) ||
-               system.surface_normal_method isa ColorfieldSurfaceNormal
+        return surface_normal_method(system) isa ColorfieldSurfaceNormal
     end
 
     if requires_color_check
+        fluid_ids = findall(contributes_to_colorfield, systems)
+        boundary_ids = findall(system -> system isa
+                                         WallBoundarySystem{<:BoundaryModelDummyParticles} ||
+                                         system isa
+                                         RigidBodySystem{<:BoundaryModelDummyParticles},
+                               systems)
+        normal_fluid_ids = filter(i -> surface_normal_method(systems[i]) isa
+                                       ColorfieldSurfaceNormal, fluid_ids)
+
+        participant_ids = (fluid_ids..., boundary_ids...)
+        if all(i -> iszero(systems[i].cache.color), participant_ids)
+            throw(ArgumentError("When `ColorfieldSurfaceNormal` is used, at least one participating system must have a color different from 0."))
+        end
 
-        # Systems that contribute to the colorfield/contact logic.
-        system_ids = findall(system -> (system isa AbstractFluidSystem &&
-                                        !(system isa ParticlePackingSystem)) ||
-                                       system isa WallBoundarySystem ||
-                                       system isa
-                                       RigidBodySystem{<:BoundaryModelDummyParticles},
-                             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."))
+        if !isempty(boundary_ids) &&
+           any(i -> iszero(systems[i].cache.color), normal_fluid_ids) &&
+           all(i -> iszero(systems[i].cache.color), boundary_ids)
+            throw(ArgumentError("A fluid with `ColorfieldSurfaceNormal` and color 0 requires " *
+                                "a nonzero dummy-boundary color for contact detection."))
         end
     end
 end

src/preprocessing/particle_packing/system.jl

@@ -100,6 +100,20 @@ struct ParticlePackingSystem{S, F, NDIMS, ELTYPE <: Real, PR, C, AV,
     end
 end
 
+@inline surface_normal_method(::ParticlePackingSystem) = nothing
+@inline contributes_to_colorfield(::ParticlePackingSystem) = false
+
+function check_configuration(system::ParticlePackingSystem, systems, nhs)
+    return nothing
+end
+
+function calc_normal!(system::AbstractFluidSystem,
+                      neighbor_system::ParticlePackingSystem,
+                      u_system, v, v_neighbor_system, u_neighbor_system, semi,
+                      surface_normal_method, neighbor_surface_normal_method)
+    return system
+end
+
 function ParticlePackingSystem(shape::InitialCondition;
                                signed_distance_field::Union{SignedDistanceField, Nothing},
                                smoothing_kernel=SchoenbergQuinticSplineKernel{ndims(shape)}(),

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/boundary/wall_boundary/system.jl

@@ -14,10 +14,9 @@ The interaction between fluid and boundary particles is specified by the boundar
 - `prescribed_motion`: For moving boundaries, a [`PrescribedMotion`](@ref) can be passed.
 - `adhesion_coefficient`: Coefficient specifying the adhesion of a fluid to the surface.
    Note: currently it is assumed that all fluids have the same adhesion coefficient.
-- `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
-   reads the resulting boundary colorfield to detect wall contact.
+- `color_value`: Integer scalar used with [`BoundaryModelDummyParticles`](@ref) and
+   [`ColorfieldSurfaceNormal`](@ref). Fluid-boundary normal evaluation uses its magnitude
+   in the smoothed boundary color field that detects wall contact.
 """
 struct WallBoundarySystem{BM, ELTYPE <: Real, NDIMS, IC, CO, M, IM,
                           CA} <: AbstractBoundarySystem{NDIMS}

src/schemes/fluid/entropically_damped_sph/system.jl

@@ -50,14 +50,16 @@ 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)
-- `reference_particle_spacing`: The reference particle spacing used for weighting values at the boundary,
-                                which currently is only needed when using surface tension.
-- `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
-                                reads the resulting boundary colorfield to detect wall contact.
+- `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`: Reference spacing used by support-based normal validity
+                                checks. It is required when using a surface-normal method.
+- `color_value`:                Integer scalar used by [`ColorfieldSurfaceNormal`](@ref).
+                                Interacting fluid values define fluid-fluid color gradients;
+                                dummy-boundary contact detection also uses the fluid value.
 
 """
 struct EntropicallyDampedSPHSystem{NDIMS, ELTYPE <: Real, IC, M, DC, K, V, COR, PF, TV,
@@ -119,12 +121,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 +309,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
@@ -254,6 +257,9 @@ end
     return nothing
 end
 
+@inline contributes_to_colorfield(system) = false
+@inline contributes_to_colorfield(::AbstractFluidSystem) = true
+
 function restart_u(system::AbstractFluidSystem, data)
     inactive_coords = convert(coordinates_eltype(system), 1e16)
     coords_total = fill(inactive_coords, u_nvariables(system),
@@ -297,7 +303,7 @@ function restart_v(system::AbstractFluidSystem, data)
 end
 
 function check_configuration(fluid_system::AbstractFluidSystem, systems, nhs)
-    if !(fluid_system isa ParticlePackingSystem) && !isnothing(fluid_system.surface_tension)
+    if !isnothing(fluid_system.surface_tension)
         foreach_system(systems) do neighbor
             if neighbor isa AbstractFluidSystem &&
                isnothing(fluid_system.surface_tension) &&

src/schemes/fluid/implicit_incompressible_sph/system.jl

@@ -4,7 +4,7 @@
                                     viscosity=nothing,
                                     acceleration=ntuple(_ -> 0.0, ndims(smoothing_kernel)),
                                     omega=0.5, max_error=0.1, min_iterations=2,
-                                    max_iterations=20, time_step)
+                                     max_iterations=20, time_step, color_value=1)
 
 System for particles of a fluid.
 The system employs implicit incompressible SPH (IISPH), iteratively solving a linear system
@@ -30,6 +30,8 @@ See [Implicit Incompressible SPH](@ref iisph) for more details on the method.
 - `min_iterations = 2`:         Minimum number of iterations in the relaxed Jacobi scheme, independent from the termination condition
 - `max_iterations = 20`:        Maximum number of iterations in the relaxed Jacobi scheme, independent from the termination condition
 - `time_step`:                  Time step size used for the simulation
+- `color_value`:                Integer scalar contributed to fluid-fluid color gradients when
+                                interacting with a fluid using [`ColorfieldSurfaceNormal`](@ref).
 """
 struct ImplicitIncompressibleSPHSystem{NDIMS, ELTYPE <: Real, ARRAY1D, ARRAY2D,
                                        IC, K, V, PF, C} <: AbstractFluidSystem{NDIMS}
@@ -71,7 +73,7 @@ function ImplicitIncompressibleSPHSystem(initial_condition; smoothing_kernel,
                                                              ndims(smoothing_kernel)),
                                          omega=0.5, max_error=0.1, min_iterations=2,
                                          max_iterations=20, time_step,
-                                         artificial_sound_speed=1000.0)
+                                         artificial_sound_speed=1000.0, color_value=1)
     particle_refinement = nothing # TODO
     surface_tension = nothing # TODO
 
@@ -125,7 +127,8 @@ function ImplicitIncompressibleSPHSystem(initial_condition; smoothing_kernel,
 
     cache = (;
              create_cache_refinement(initial_condition, particle_refinement,
-                                     smoothing_length)...,)
+                                     smoothing_length)...,
+             color=Int(color_value))
 
     return ImplicitIncompressibleSPHSystem(initial_condition, mass, pressure,
                                            smoothing_kernel, smoothing_length,