Use conservative Akinci pair forces across resolutions

NEWS.md

@@ -4,6 +4,36 @@ 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
+
+### API Changes
+
+- `SurfaceTensionAkinci` now requires `AkinciFreeSurfaceCorrection`, which implements the
+  neighborhood-deficiency factor required by the published combined force.
+- Akinci cohesion and adhesion kernels now use dimensionally consistent, integral-matched
+  normalizations in 2D. To preserve previous pairwise kernel contributions at compact-support
+  radius `h_c`, multiply the surface-tension coefficient by `627 / (790 * h_c)` and the
+  adhesion coefficient by `42 / (65 * h_c)`; migrated coefficients are resolution-independent.
+
+### Features
+
+- Added `AkinciFreeSurfaceCorrection` support to `EntropicallyDampedSPHSystem` for both
+  continuity and summation density.
+- Akinci fluid-fluid forces now use a symmetric pair formulation for differing particle masses
+  and smoothing lengths while preserving the equal-resolution result.
+
+### 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.
+- Fixed `AkinciFreeSurfaceCorrection` with `ContinuityDensity` by reconstructing the
+  neighborhood-deficiency density used by the correction and the `SurfaceTensionAkinci`
+  color-field normals without changing the integrated pressure density. Dummy boundary masses
+  are included so wall particles are not misclassified as a free surface.
+
 ## Version 0.5.3
 
 ### Features

docs/src/systems/entropically_damped_sph.md

@@ -35,6 +35,11 @@ The viscosity parameter ``\eta_a`` for a particle ``a`` is given as
 where it is found in the numerical experiments of [Ramachandran (2019)](@cite Ramachandran2019) that ``\alpha = 0.5``
 is a good choice for a wide range of Reynolds numbers (0.0125 to 10000).
 
+[`AkinciFreeSurfaceCorrection`](@ref) can be used with either density calculator. It scales
+viscosity and Akinci surface-tension forces without modifying EDAC pressure acceleration or
+pressure evolution. With [`ContinuityDensity`](@ref), an auxiliary summation density is used to
+measure free-surface neighborhood deficiency.
+
 !!! note
     > The EDAC formulation keeps the density constant and this eliminates the need for the continuity equation
     > or the use of a summation density to find the pressure. However, in SPH discretizations, ``m/\rho``

examples/fluid/falling_water_spheres_2d.jl

@@ -61,13 +61,16 @@ viscosity = ArtificialViscosityMonaghan(; alpha, beta=0.0)
 density_diffusion = DensityDiffusionAntuono(delta=0.1)
 surface_tension_coefficient = 0.05
 surface_tension = SurfaceTensionAkinci(; surface_tension_coefficient)
+surface_tension_correction = surface_tension isa SurfaceTensionAkinci ?
+                             AkinciFreeSurfaceCorrection(fluid_density) : nothing
 
 sphere_surface_tension = EntropicallyDampedSPHSystem(sphere1;
                                                      smoothing_kernel=fluid_smoothing_kernel,
                                                      smoothing_length=fluid_smoothing_length,
                                                      sound_speed, viscosity,
                                                      density_calculator=ContinuityDensity(),
                                                      acceleration, surface_tension,
+                                                     correction=surface_tension_correction,
                                                      reference_particle_spacing=fluid_particle_spacing)
 
 sphere = WeaklyCompressibleSPHSystem(sphere2; smoothing_kernel=fluid_smoothing_kernel,

src/general/corrections.jl

@@ -3,12 +3,12 @@
     AkinciFreeSurfaceCorrection(rho0)
 
 Free surface correction according to [Akinci et al. (2013)](@cite Akinci2013).
-At a free surface, the mean density is typically lower than the reference density,
-resulting in reduced surface tension and viscosity forces.
-The free surface correction adjusts the viscosity, pressure, and surface tension forces
-near free surfaces to counter this effect.
-It's important to note that this correlation is unphysical and serves as an approximation.
-The computation time added by this method is about 2--3%.
+At a free surface, particle-neighborhood deficiency reduces the pairwise surface tension and
+viscosity forces. This correction applies the symmetrized factor from Equation 4 to the combined
+cohesion and curvature force in Equation 5 and, as specified in Section 4, to viscosity. It does
+not modify pressure forces.
+The published implementation reports about 2--3% overhead for evaluating the correction from
+an already available density estimate.
 
 Mathematically the idea is quite simple. If we have an SPH particle in the middle of a volume
 at rest, its density will be identical to the rest density ``\rho_0``. If we now consider an SPH
@@ -18,7 +18,15 @@ the surface, which will result in a lower density. If we calculate the correctio
 k = \rho_0/\rho_\text{mean},
 ```
 this value will be about ~1.5 for particles at the free surface and can then be used to increase
-the pressure and viscosity accordingly.
+the surface tension and viscosity forces accordingly.
+
+With [`SummationDensity`](@ref), the correction uses the current density directly. When used with
+[`ContinuityDensity`](@ref) in a [`WeaklyCompressibleSPHSystem`](@ref) or
+[`EntropicallyDampedSPHSystem`](@ref), TrixiParticles.jl reconstructs an auxiliary summation
+density for this correction and for the color-field normals of [`SurfaceTensionAkinci`](@ref),
+while pressure continues to use the integrated density. This preserves the
+particle-neighborhood-deficiency behavior of the published model at the cost of one additional
+density summation per update stage.
 
 # Arguments
 - `rho0`: Rest density.
@@ -40,7 +48,7 @@ end
     rho_mean = (rho_a + rho_b) / 2
     k = correction.rho0 / rho_mean
 
-    # Viscosity, pressure, surface_tension
+    # Equation 5 applies `k` to surface tension; Section 4 also applies it to viscosity.
     return k, 1, k
 end
 
@@ -458,6 +466,10 @@ end
 
 create_cache_correction(correction, density, NDIMS, nparticles) = (;)
 
+function create_cache_correction(::AkinciFreeSurfaceCorrection, density, NDIMS, n_particles)
+    return (; kernel_summation_density=similar(density))
+end
+
 function create_cache_correction(::ShepardKernelCorrection, density, NDIMS, n_particles)
     return (; kernel_correction_coefficient=similar(density))
 end

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/rhs.jl

@@ -42,6 +42,14 @@ function interact!(dv, v_particle_system, u_particle_system,
         rho_a = @inbounds current_density(v_particle_system, particle_system, particle)
         rho_b = @inbounds current_density(v_neighbor_system, neighbor_system, neighbor)
 
+        correction_rho_a = correction_density(correction, particle_system, particle, rho_a)
+        correction_rho_b = correction_density(correction, neighbor_system, neighbor, rho_b)
+        (viscosity_correction, pressure_correction,
+         surface_tension_correction) = free_surface_correction(correction,
+                                                               particle_system,
+                                                               correction_rho_a,
+                                                               correction_rho_b)
+
         v_a = @inbounds current_velocity(v_particle_system, particle_system, particle)
         v_b = @inbounds current_velocity(v_neighbor_system, neighbor_system, neighbor)
 
@@ -65,12 +73,12 @@ function interact!(dv, v_particle_system, u_particle_system,
                                             rho_b, pos_diff, distance, grad_kernel,
                                             correction)
 
-        dv_particle = Ref(dv_pressure)
+        dv_particle = Ref(pressure_correction * dv_pressure)
         @inbounds dv_viscosity!(dv_particle, particle_system, neighbor_system,
                                 v_particle_system, v_neighbor_system,
                                 particle, neighbor, pos_diff, distance,
                                 sound_speed, m_a, m_b, rho_a, rho_b,
-                                v_a, v_b, grad_kernel)
+                                v_a, v_b, grad_kernel, viscosity_correction)
 
         # Extra terms in the momentum equation when using a shifting technique
         @inbounds dv_shifting!(dv_particle, shifting_technique(particle_system),
@@ -83,7 +91,8 @@ function interact!(dv, v_particle_system, u_particle_system,
                                          surface_tension_b,
                                          particle_system, neighbor_system,
                                          particle, neighbor, pos_diff, distance,
-                                         rho_a, rho_b, grad_kernel, 1)
+                                         rho_a, rho_b, grad_kernel,
+                                         surface_tension_correction)
 
         @inbounds adhesion_force!(dv_particle, surface_tension_a, particle_system,
                                   neighbor_system,