Generalize free-surface detection

NEWS.md

@@ -4,9 +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.3
-
-### Features
+## Version 0.6.0
+
+### API Changes
+
+- Renamed the fluid-system keyword `surface_normal_method` to `surface_method` and added
+  detection-only surface methods. The old constructor keyword and accessor are deprecated.
+
+## Version 0.5.3
+
+### Features
 
 - Added normal vectors to `InitialCondition`, with automatic computation for boundary
   particles generated by `RectangularTank` and `SphereShape` (#1036).

docs/src/systems/fluid.md

@@ -211,7 +211,7 @@ Pages = [joinpath("general", "corrections.jl")]
 
 ---
 
-## [Surface Normals](@id surface_normal)
+## [Surface Detection And Normals](@id surface_normal)
 
 ### Overview of surface normal calculation in SPH
 
@@ -248,6 +248,28 @@ The calculated normals are normalized to unit vectors:
 
 Normalization ensures that the magnitude of the normals does not bias the curvature calculations or the resulting surface tension forces.
 
+#### Surface methods and activity
+
+Fluid systems configure interface geometry with the `surface_method` keyword. Every surface
+method computes a smooth `surface_activity`. Methods derived from
+`AbstractSurfaceNormalMethod` additionally provide a normal, so normal calculation always
+includes detection.
+
+`ColorfieldSurfaceDetection` computes activity only. `ColorfieldSurfaceNormal` uses the same
+colorfield accumulation and additionally filters and stores the gradient as a surface normal.
+Different `color_value`s detect interfaces between represented liquids. A constant nonzero
+color detects a free surface because its kernel support ends at the unrepresented exterior.
+Equal colors do not create an internal interface.
+
+`surface_activity` is available in particle VTK output and as a custom quantity for
+`SolutionSavingCallback` and `PostprocessCallback`. `surf_normal` is written only for
+normal-capable methods.
+
+Point and plane interpolation evaluate the same colorfield gradient. With `cut_off_bnd=true`,
+kernel-weighted color contributions also determine whether a point belongs to the reference
+phase. This prevents extrapolation through both free surfaces and interfaces with another
+liquid while retaining the existing solid-boundary cutoff.
+
 #### Handling noise and errors in normal calculation
 
 In regions distant from the interface, the calculated normals may be small or inaccurate due to the

examples/fluid/sphere_surface_tension_2d.jl

@@ -53,7 +53,7 @@ fluid_system = EntropicallyDampedSPHSystem(fluid; smoothing_kernel=fluid_smoothi
                                            density_calculator=ContinuityDensity(),
                                            reference_particle_spacing=particle_spacing,
                                            acceleration=zeros(length(fluid_size)),
-                                           surface_normal_method=ColorfieldSurfaceNormal(),
+                                           surface_method=ColorfieldSurfaceNormal(),
                                            surface_tension=SurfaceTensionMorris(surface_tension_coefficient=50 *
                                                                                                             0.0728))
 

src/TrixiParticles.jl

@@ -115,7 +115,9 @@ export interpolate_line, interpolate_points, interpolate_plane_3d, interpolate_p
        interpolate_plane_2d_vtk
 export SurfaceTensionAkinci, CohesionForceAkinci, SurfaceTensionMorris,
        SurfaceTensionMomentumMorris
-export ColorfieldSurfaceNormal
+export AbstractSurfaceMethod, AbstractSurfaceNormalMethod, ColorfieldSurfaceDetection,
+       ColorfieldSurfaceNormal, surface_method, computes_surface_normal, surface_activity,
+       surface_normal
 export SymplecticPositionVerlet
 export coordinates_eltype
 

src/general/custom_quantities.jl

@@ -141,3 +141,24 @@ end
 function avg_density(system, dv_ode, du_ode, v_ode, u_ode, semi, t)
     return NaN
 end
+
+"""
+    surface_activity
+
+Return the per-particle smooth surface activity, or `nothing` for systems without a
+configured surface method.
+"""
+function surface_activity(system::AbstractFluidSystem, dv_ode, du_ode, v_ode, u_ode,
+                          semi, t)
+    hasproperty(system.cache, :surface_activity) || return nothing
+    return view(system.cache.surface_activity, each_active_particle(system))
+end
+
+surface_activity(system, dv_ode, du_ode, v_ode, u_ode, semi, t) = nothing
+
+function surface_normal(system::AbstractFluidSystem, dv_ode, du_ode, v_ode, u_ode, semi, t)
+    computes_surface_normal(surface_method(system)) || return nothing
+    return view(system.cache.surface_normal, :, each_active_particle(system))
+end
+
+surface_normal(system, dv_ode, du_ode, v_ode, u_ode, semi, t) = nothing