From b0aef6cad859244be308225b76554935b6ad3919 Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Thu, 6 Aug 2026 11:40:03 +0200 Subject: [PATCH 01/13] Harden surface tension model configuration --- NEWS.md | 8 ++ docs/src/systems/fluid.md | 23 +++- .../boundary/wall_boundary/dummy_particles.jl | 2 +- .../fluid/entropically_damped_sph/system.jl | 12 +- src/schemes/fluid/fluid.jl | 11 +- src/schemes/fluid/surface_normal_sph.jl | 13 +- src/schemes/fluid/surface_tension.jl | 77 +++++++---- .../fluid/weakly_compressible_sph/system.jl | 12 +- test/schemes/fluid/surface_tension.jl | 126 +++++++++++++++++- 9 files changed, 236 insertions(+), 48 deletions(-) diff --git a/NEWS.md b/NEWS.md index ad3e44c715..7168811c57 100644 --- a/NEWS.md +++ b/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, avoiding unnecessary + surface-normal allocation for `CohesionForceAkinci`, and stabilizing Akinci adhesion evaluation + near the compact-support boundary. + ## Version 0.5.3 ### Features diff --git a/docs/src/systems/fluid.md b/docs/src/systems/fluid.md index c9b9860ab9..63bccfc66c 100644 --- a/docs/src/systems/fluid.md +++ b/docs/src/systems/fluid.md @@ -215,9 +215,9 @@ 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 provide the directionality of forces acting at the fluid interface. They are +used by the full Akinci model and both Morris models, but not by the cohesion-only Akinci model. +They are calculated based on the particle properties and their spatial distribution. #### Color field and gradient-based surface normals @@ -289,6 +289,23 @@ 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. + +!!! warning "Akinci coefficients in two dimensions" + The Akinci cohesion and adhesion kernels use the normalization published for the + three-dimensional model. In two-dimensional simulations, their coefficients are empirical + numerical parameters rather than resolution-independent physical values in N/m. Recheck the + coefficient when changing particle spacing or smoothing length. + ### [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: diff --git a/src/schemes/boundary/wall_boundary/dummy_particles.jl b/src/schemes/boundary/wall_boundary/dummy_particles.jl index fcde467f2e..437610d9a8 100644 --- a/src/schemes/boundary/wall_boundary/dummy_particles.jl +++ b/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 diff --git a/src/schemes/fluid/entropically_damped_sph/system.jl b/src/schemes/fluid/entropically_damped_sph/system.jl index 6e03cf7837..1ea4b5b74e 100644 --- a/src/schemes/fluid/entropically_damped_sph/system.jl +++ b/src/schemes/fluid/entropically_damped_sph/system.jl @@ -51,9 +51,10 @@ See [Entropically Damped Artificial Compressibility for SPH](@ref edac) for more 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) + (default: no surface normal method or `ColorfieldSurfaceNormal()` + if 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 +120,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 && diff --git a/src/schemes/fluid/fluid.jl b/src/schemes/fluid/fluid.jl index fb49a65359..d94943f7bc 100644 --- a/src/schemes/fluid/fluid.jl +++ b/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 diff --git a/src/schemes/fluid/surface_normal_sph.jl b/src/schemes/fluid/surface_normal_sph.jl index adbc9d7dbe..3562b93ddf 100644 --- a/src/schemes/fluid/surface_normal_sph.jl +++ b/src/schemes/fluid/surface_normal_sph.jl @@ -17,8 +17,17 @@ 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 = promote(boundary_contact_threshold, interface_threshold, + ideal_density_threshold) + 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) diff --git a/src/schemes/fluid/surface_tension.jl b/src/schemes/fluid/surface_tension.jl index 5656e95e12..5673125016 100644 --- a/src/schemes/fluid/surface_tension.jl +++ b/src/schemes/fluid/surface_tension.jl @@ -1,22 +1,38 @@ abstract type AbstractSurfaceTension end abstract type AkinciTypeSurfaceTension <: AbstractSurfaceTension end +function validate_surface_tension_coefficient(surface_tension_coefficient) + if !(surface_tension_coefficient isa Real) || + !isfinite(surface_tension_coefficient) || surface_tension_coefficient < 0 + throw(ArgumentError("`surface_tension_coefficient` must be a finite, non-negative real number")) + end + + return surface_tension_coefficient +end + @doc raw""" CohesionForceAkinci(surface_tension_coefficient=1.0) This model only implements the cohesion force of the Akinci [Akinci2013](@cite) surface tension model. +It does not require a surface-normal method. + +The published Akinci cohesion kernel uses a three-dimensional normalization. In two-dimensional +simulations, `surface_tension_coefficient` is therefore an empirical numerical parameter and +may need to be adjusted when changing the resolution. See [`surface_tension`](@ref) for more details. # Keywords -- `surface_tension_coefficient=1.0`: Modifies the intensity of the surface tension-induced force, - enabling the tuning of the fluid's surface tension properties within the simulation. +- `surface_tension_coefficient=1.0`: Finite, non-negative coefficient modifying the + fluid-fluid cohesion force. Zero disables this force; wall adhesion is controlled by the + boundary's `adhesion_coefficient`. """ -struct CohesionForceAkinci{ELTYPE} <: AkinciTypeSurfaceTension +struct CohesionForceAkinci{ELTYPE <: Real} <: AkinciTypeSurfaceTension surface_tension_coefficient::ELTYPE function CohesionForceAkinci(; surface_tension_coefficient=1.0) - new{typeof(surface_tension_coefficient)}(surface_tension_coefficient) + coefficient = validate_surface_tension_coefficient(surface_tension_coefficient) + new{typeof(coefficient)}(coefficient) end end @@ -28,18 +44,22 @@ principles outlined by Akinci [Akinci2013](@cite). This model is instrumental in behaviors of fluid surfaces, such as droplet formation and the dynamics of merging or separation, by utilizing intra-particle forces. +The published Akinci cohesion kernel uses a three-dimensional normalization. In two-dimensional +simulations, `surface_tension_coefficient` is therefore an empirical numerical parameter and +may need to be adjusted when changing the resolution. + See [`surface_tension`](@ref) for more details. # Keywords -- `surface_tension_coefficient=1.0`: A parameter to adjust the magnitude of - surface tension forces, facilitating the fine-tuning of how surface tension phenomena - are represented in the simulation. +- `surface_tension_coefficient=1.0`: Finite, non-negative coefficient adjusting the + magnitude of surface tension forces. Zero disables the fluid-fluid force. """ -struct SurfaceTensionAkinci{ELTYPE} <: AkinciTypeSurfaceTension +struct SurfaceTensionAkinci{ELTYPE <: Real} <: AkinciTypeSurfaceTension surface_tension_coefficient::ELTYPE function SurfaceTensionAkinci(; surface_tension_coefficient=1.0) - new{typeof(surface_tension_coefficient)}(surface_tension_coefficient) + coefficient = validate_surface_tension_coefficient(surface_tension_coefficient) + new{typeof(coefficient)}(coefficient) end end @@ -55,14 +75,15 @@ See [`surface_tension`](@ref) for more details. # Keywords -- `surface_tension_coefficient=1.0`: Adjusts the magnitude of the surface tension - forces, enabling tuning of fluid surface behaviors in simulations. +- `surface_tension_coefficient=1.0`: Finite, non-negative coefficient adjusting the + magnitude of surface tension forces. Zero disables the force. """ -struct SurfaceTensionMorris{ELTYPE} <: AbstractSurfaceTension +struct SurfaceTensionMorris{ELTYPE <: Real} <: AbstractSurfaceTension surface_tension_coefficient::ELTYPE function SurfaceTensionMorris(; surface_tension_coefficient=1.0) - new{typeof(surface_tension_coefficient)}(surface_tension_coefficient) + coefficient = validate_surface_tension_coefficient(surface_tension_coefficient) + new{typeof(coefficient)}(coefficient) end end @@ -87,17 +108,24 @@ numerical adjustments at higher resolutions. See [`surface_tension`](@ref) for more details. # Keywords -- `surface_tension_coefficient=1.0`: A parameter to adjust the strength of surface tension - forces, allowing fine-tuning to replicate physical behavior. +- `surface_tension_coefficient=1.0`: Finite, non-negative coefficient adjusting the + strength of surface tension forces. Zero disables the force. """ -struct SurfaceTensionMomentumMorris{ELTYPE} <: AbstractSurfaceTension +struct SurfaceTensionMomentumMorris{ELTYPE <: Real} <: AbstractSurfaceTension surface_tension_coefficient::ELTYPE function SurfaceTensionMomentumMorris(; surface_tension_coefficient=1.0) - new{typeof(surface_tension_coefficient)}(surface_tension_coefficient) + coefficient = validate_surface_tension_coefficient(surface_tension_coefficient) + new{typeof(coefficient)}(coefficient) end end +# Surface-model capabilities are expressed through dispatch so that constructors do not need +# to duplicate concrete model checks. +@inline requires_surface_normal(::Nothing) = false +@inline requires_surface_normal(::CohesionForceAkinci) = false +@inline requires_surface_normal(::Any) = true + function create_cache_surface_tension(::SurfaceTensionMomentumMorris, ELTYPE, NDIMS, nparticles) delta_s = Array{ELTYPE, 1}(undef, nparticles) @@ -141,18 +169,17 @@ end @inline function adhesion_force_akinci(surface_tension, support_radius, m_b, pos_diff, distance, adhesion_coefficient) - - # The neighborhood search has an `<=` check, but for `distance == support_radius` - # the term inside the parentheses might be very slightly negative, causing an error with `^0.25`. - # TODO Change this in the neighborhood search? - # See https://github.com/trixi-framework/PointNeighbors.jl/issues/19 distance >= support_radius && return zero(pos_diff) distance <= 0.5 * support_radius && return zero(pos_diff) - # Eq. 7 - A = 0.007 / support_radius^3.25 * - (-4 * distance^2 / support_radius + 6 * distance - 2 * support_radius)^0.25 + # Eq. 7. The factored radicand avoids cancellation close to the support boundary. + radicand = 2 * (2 * distance - support_radius) * + (support_radius - distance) / support_radius + fourth_root = sqrt(sqrt(max(zero(radicand), radicand))) + normalization = convert(typeof(support_radius), 0.007) / + (support_radius^3 * sqrt(sqrt(support_radius))) + A = normalization * fourth_root # Eq. 6 in acceleration form with `m_b` being the boundary mass calculated as # `m_b = rho_0 * volume` (Akinci boundary condition treatment) diff --git a/src/schemes/fluid/weakly_compressible_sph/system.jl b/src/schemes/fluid/weakly_compressible_sph/system.jl index eea0607d7d..c97c829824 100644 --- a/src/schemes/fluid/weakly_compressible_sph/system.jl +++ b/src/schemes/fluid/weakly_compressible_sph/system.jl @@ -54,9 +54,10 @@ See [Weakly Compressible SPH](@ref wcsph) for more details on the method. 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) + (default: no surface normal method or `ColorfieldSurfaceNormal()` + if 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 @@ -130,12 +131,11 @@ function WeaklyCompressibleSPHSystem(initial_condition; smoothing_kernel, throw(ArgumentError("`ShepardKernelCorrection` cannot be used with `ContinuityDensity`")) 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 pressure_acceleration = choose_pressure_acceleration_formulation(pressure_acceleration, diff --git a/test/schemes/fluid/surface_tension.jl b/test/schemes/fluid/surface_tension.jl index 7fe8abbd97..d470da0ba9 100644 --- a/test/schemes/fluid/surface_tension.jl +++ b/test/schemes/fluid/surface_tension.jl @@ -1,5 +1,119 @@ - @testset verbose=true "Surface Tension" begin + @testset "constructors and capabilities" begin + constructors = (CohesionForceAkinci, SurfaceTensionAkinci, + SurfaceTensionMorris, SurfaceTensionMomentumMorris) + + for constructor in constructors + model = constructor(surface_tension_coefficient=0.5f0) + @test model.surface_tension_coefficient === 0.5f0 + @test iszero(constructor(surface_tension_coefficient=0).surface_tension_coefficient) + + for coefficient in (-1.0, NaN, Inf, -Inf, 1.0im, "invalid") + @test_throws ArgumentError constructor(surface_tension_coefficient=coefficient) + end + end + + @test !TrixiParticles.requires_surface_normal(nothing) + @test !TrixiParticles.requires_surface_normal(CohesionForceAkinci()) + @test TrixiParticles.requires_surface_normal(SurfaceTensionAkinci()) + @test TrixiParticles.requires_surface_normal(SurfaceTensionMorris()) + @test TrixiParticles.requires_surface_normal(SurfaceTensionMomentumMorris()) + + normal_method = ColorfieldSurfaceNormal(boundary_contact_threshold=1, + interface_threshold=0.1f0, + ideal_density_threshold=0.25) + @test normal_method isa ColorfieldSurfaceNormal{Float64} + @test ColorfieldSurfaceNormal(boundary_contact_threshold=0.1f0, + interface_threshold=0.01f0, + ideal_density_threshold=0.0f0) isa + ColorfieldSurfaceNormal{Float32} + end + + @testset "cohesion-only systems do not require normals" begin + coordinates = [0.0 1.0; + 0.0 0.0] + initial_condition = InitialCondition(; coordinates, density=ones(2), + particle_spacing=1.0) + smoothing_kernel = WendlandC2Kernel{2}() + smoothing_length = 1.0 + surface_tension = CohesionForceAkinci(surface_tension_coefficient=0.1) + + wcsph = WeaklyCompressibleSPHSystem(initial_condition; smoothing_kernel, + smoothing_length, + density_calculator=SummationDensity(), + state_equation=StateEquationCole(sound_speed=10.0, + reference_density=1.0, + exponent=1), + surface_tension) + edac = EntropicallyDampedSPHSystem(initial_condition; smoothing_kernel, + smoothing_length, sound_speed=10.0, + density_calculator=SummationDensity(), + surface_tension) + + for system in (wcsph, edac) + @test isnothing(system.surface_normal_method) + @test !haskey(system.cache, :surface_normal) + @test !haskey(system.cache, :neighbor_count) + @test !haskey(system.cache, :reference_particle_spacing) + end + + @test_throws ArgumentError WeaklyCompressibleSPHSystem(initial_condition; + smoothing_kernel, + smoothing_length, + density_calculator=SummationDensity(), + state_equation=StateEquationCole(sound_speed=10.0, + reference_density=1.0, + exponent=1), + surface_tension=SurfaceTensionAkinci()) + @test_throws ArgumentError EntropicallyDampedSPHSystem(initial_condition; + smoothing_kernel, + smoothing_length, + sound_speed=10.0, + density_calculator=SummationDensity(), + surface_tension=SurfaceTensionAkinci()) + + full_akinci = WeaklyCompressibleSPHSystem(initial_condition; smoothing_kernel, + smoothing_length, + density_calculator=SummationDensity(), + state_equation=StateEquationCole(sound_speed=10.0, + reference_density=1.0, + exponent=1), + surface_tension=SurfaceTensionAkinci(), + reference_particle_spacing=1.0) + @test full_akinci.surface_normal_method isa ColorfieldSurfaceNormal + @test haskey(full_akinci.cache, :surface_normal) + end + + @testset "zero Morris coefficient does not restrict the time step" begin + function calculate_initial_dt(surface_tension) + initial_condition = InitialCondition(; coordinates=[0.0 1.0; 0.0 0.0], + density=ones(2), particle_spacing=1.0) + reference_particle_spacing = isnothing(surface_tension) ? 0 : 1.0 + system = WeaklyCompressibleSPHSystem(initial_condition; + smoothing_kernel=WendlandC2Kernel{2}(), + smoothing_length=1.0, + density_calculator=SummationDensity(), + state_equation=StateEquationCole(sound_speed=10.0, + reference_density=1.0, + exponent=1), + surface_tension, + reference_particle_spacing) + semi = Semidiscretization(system) + ode = semidiscretize(semi, (0.0, 0.1)) + v_ode, u_ode = ode.u0.x + return TrixiParticles.calculate_dt(v_ode, u_ode, 0.25, semi.systems[1], semi) + end + + dt_without_surface_tension = calculate_initial_dt(nothing) + dt_with_zero_csf = calculate_initial_dt(SurfaceTensionMorris(; + surface_tension_coefficient=0.0)) + dt_with_zero_css = calculate_initial_dt(SurfaceTensionMomentumMorris(; + surface_tension_coefficient=0.0)) + + @test dt_with_zero_csf == dt_without_surface_tension + @test dt_with_zero_css == dt_without_surface_tension + end + @testset verbose=true "`cohesion_force_akinci`" begin surface_tension = SurfaceTensionAkinci(surface_tension_coefficient=1.0) support_radius = 1.0 @@ -88,6 +202,16 @@ test_distance @test isapprox(zero[1], 0.0, atol=6e-15) @test isapprox(zero[2], 0.0, atol=6e-15) + + support_radius_f32 = 15.594092f0 + distance_f32 = prevfloat(support_radius_f32) + near_support = TrixiParticles.adhesion_force_akinci(surface_tension, + support_radius_f32, 1.0f0, + Float32[1, 0], distance_f32, + 1.0f0) + @test eltype(near_support) == Float32 + @test all(isfinite, near_support) + @test 0 < norm(near_support) < eps(Float32) end @testset "compute_stress_tensors! (MomentumMorris)" begin From fd596fcb1fcc4f2f0c6e24694a6cd5084eaedcaa Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Fri, 14 Aug 2026 14:33:37 +0200 Subject: [PATCH 02/13] Remove redundant surface tension type bounds --- src/schemes/fluid/surface_tension.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/schemes/fluid/surface_tension.jl b/src/schemes/fluid/surface_tension.jl index 5673125016..53e6151363 100644 --- a/src/schemes/fluid/surface_tension.jl +++ b/src/schemes/fluid/surface_tension.jl @@ -27,7 +27,7 @@ See [`surface_tension`](@ref) for more details. fluid-fluid cohesion force. Zero disables this force; wall adhesion is controlled by the boundary's `adhesion_coefficient`. """ -struct CohesionForceAkinci{ELTYPE <: Real} <: AkinciTypeSurfaceTension +struct CohesionForceAkinci{ELTYPE} <: AkinciTypeSurfaceTension surface_tension_coefficient::ELTYPE function CohesionForceAkinci(; surface_tension_coefficient=1.0) @@ -54,7 +54,7 @@ See [`surface_tension`](@ref) for more details. - `surface_tension_coefficient=1.0`: Finite, non-negative coefficient adjusting the magnitude of surface tension forces. Zero disables the fluid-fluid force. """ -struct SurfaceTensionAkinci{ELTYPE <: Real} <: AkinciTypeSurfaceTension +struct SurfaceTensionAkinci{ELTYPE} <: AkinciTypeSurfaceTension surface_tension_coefficient::ELTYPE function SurfaceTensionAkinci(; surface_tension_coefficient=1.0) @@ -78,7 +78,7 @@ See [`surface_tension`](@ref) for more details. - `surface_tension_coefficient=1.0`: Finite, non-negative coefficient adjusting the magnitude of surface tension forces. Zero disables the force. """ -struct SurfaceTensionMorris{ELTYPE <: Real} <: AbstractSurfaceTension +struct SurfaceTensionMorris{ELTYPE} <: AbstractSurfaceTension surface_tension_coefficient::ELTYPE function SurfaceTensionMorris(; surface_tension_coefficient=1.0) @@ -111,7 +111,7 @@ See [`surface_tension`](@ref) for more details. - `surface_tension_coefficient=1.0`: Finite, non-negative coefficient adjusting the strength of surface tension forces. Zero disables the force. """ -struct SurfaceTensionMomentumMorris{ELTYPE <: Real} <: AbstractSurfaceTension +struct SurfaceTensionMomentumMorris{ELTYPE} <: AbstractSurfaceTension surface_tension_coefficient::ELTYPE function SurfaceTensionMomentumMorris(; surface_tension_coefficient=1.0) From a79806e928c5dca7c4e4b137a30ef7d60583ba42 Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Fri, 14 Aug 2026 14:51:35 +0200 Subject: [PATCH 03/13] Remove transient 2D Akinci coefficient warning --- docs/src/systems/fluid.md | 6 ------ 1 file changed, 6 deletions(-) diff --git a/docs/src/systems/fluid.md b/docs/src/systems/fluid.md index 63bccfc66c..d1f5adcc5b 100644 --- a/docs/src/systems/fluid.md +++ b/docs/src/systems/fluid.md @@ -300,12 +300,6 @@ 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. -!!! warning "Akinci coefficients in two dimensions" - The Akinci cohesion and adhesion kernels use the normalization published for the - three-dimensional model. In two-dimensional simulations, their coefficients are empirical - numerical parameters rather than resolution-independent physical values in N/m. Recheck the - coefficient when changing particle spacing or smoothing length. - ### [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: From 2134d2f1797ce78f93ec93c863756251f093d7c2 Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Fri, 14 Aug 2026 16:30:45 +0200 Subject: [PATCH 04/13] review --- NEWS.md | 7 +- docs/src/systems/fluid.md | 45 ++++- examples/fluid/falling_water_spheres_2d.jl | 3 + src/general/corrections.jl | 4 + .../fluid/entropically_damped_sph/rhs.jl | 15 +- .../fluid/entropically_damped_sph/system.jl | 23 ++- src/schemes/fluid/surface_normal_sph.jl | 13 +- src/schemes/fluid/surface_tension.jl | 95 +++++++--- .../fluid/weakly_compressible_sph/rhs.jl | 21 ++- .../fluid/weakly_compressible_sph/system.jl | 35 ++++ test/schemes/fluid/surface_normal_sph.jl | 4 + test/schemes/fluid/surface_tension.jl | 165 +++++++++++++++++- test/systems/rigid_system.jl | 1 + 13 files changed, 387 insertions(+), 44 deletions(-) diff --git a/NEWS.md b/NEWS.md index 7168811c57..49bc89f880 100644 --- a/NEWS.md +++ b/NEWS.md @@ -8,9 +8,10 @@ used in the Julia ecosystem. Notable changes will be documented in this file for ### Important Bugfixes -- Hardened surface tension model configuration by validating coefficients, avoiding unnecessary - surface-normal allocation for `CohesionForceAkinci`, and stabilizing Akinci adhesion evaluation - near the compact-support boundary. +- Hardened surface tension model configuration by validating coefficients, avoiding unnecessary + surface-normal allocation for `CohesionForceAkinci`, requiring the published free-surface + correction for the full Akinci model, stabilizing Akinci kernels across floating-point scales, + and conserving pair momentum across unequal particle masses and smoothing lengths. ## Version 0.5.3 diff --git a/docs/src/systems/fluid.md b/docs/src/systems/fluid.md index d1f5adcc5b..b518f2de85 100644 --- a/docs/src/systems/fluid.md +++ b/docs/src/systems/fluid.md @@ -240,13 +240,15 @@ where: #### Normalization of surface normals -The calculated normals are normalized to unit vectors: +Models that require only the normal direction normalize the calculated gradient: ```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. +The Morris models use this unit normal. The Akinci surface-area term instead retains the raw +gradient magnitude and multiplies it by the particle smoothing length, as required by the +published formulation. #### Handling noise and errors in normal calculation @@ -299,6 +301,8 @@ the fluid-fluid surface force. Wall adhesion is controlled independently by the 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. +The full Akinci model also requires `AkinciFreeSurfaceCorrection`, which supplies the +particle-neighborhood correction in Equations 4 and 5 of Akinci et al. ### [Akinci-based intra-particle force surface tension and wall adhesion model](@id akinci_ipf) @@ -333,14 +337,39 @@ C(r)=\frac{32}{\pi h_c^9} #### Surface area minimization force -The surface area minimization force models the curvature reduction effects, aligning particle motion to reduce the interface's total area. -It acts based on the difference in surface normals: +The surface area minimization force models curvature reduction. Let ``n_i`` denote the raw color +gradient and ``N_i=h_i n_i`` its dimensionless Akinci normal. For equal particle masses, the +published pair force is ```math -F_{\text{curvature}} = -\sigma (n_a - n_b), +F_{a\leftarrow b}^{\text{curvature}} = -\sigma m_a(N_a-N_b), ``` -where ``n_a`` and ``n_b`` are the surface normals of the interacting particles. +For equal particle masses and smoothing lengths, this is the published Akinci discretization. +For adaptive resolution, TrixiParticles forms the dimensionless normals +``N_i=h_i n_i`` and uses the symmetric pair force + +```math +F_{a\leftarrow b}^{\text{curvature}} = +-\sigma\frac{2m_am_b}{m_a+m_b}(N_a-N_b). +``` + +For a shared correction factor, the corresponding accelerations are equal and opposite after +multiplication by particle mass, and reduce to the published expression when ``m_a=m_b``. The +cohesion and curvature forces use the smaller of the two compact-support radii. This guarantees +that both directed neighborhood searches contain every active pair. Coefficients from two +different materials are averaged only as a conservative numerical default; this is not a physical +multiphase interfacial-tension law. A future multiphase model must represent the independent pair +coefficient and the pair reference density explicitly. + +Both Akinci contributions are multiplied by the required neighborhood-deficiency correction + +```math +K_{ab}=\frac{2\rho_0}{\rho_a+\rho_b}, +``` + +implemented by `AkinciFreeSurfaceCorrection`. With `ContinuityDensity`, an auxiliary summation +density is reconstructed for this correction so that missing neighbors remain observable. #### Wall adhesion force @@ -361,6 +390,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) diff --git a/examples/fluid/falling_water_spheres_2d.jl b/examples/fluid/falling_water_spheres_2d.jl index 1f0016ac6d..208c1b71bc 100644 --- a/examples/fluid/falling_water_spheres_2d.jl +++ b/examples/fluid/falling_water_spheres_2d.jl @@ -61,6 +61,8 @@ 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, @@ -68,6 +70,7 @@ sphere_surface_tension = EntropicallyDampedSPHSystem(sphere1; 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, diff --git a/src/general/corrections.jl b/src/general/corrections.jl index d97eeaf341..f6f8fcdf9d 100644 --- a/src/general/corrections.jl +++ b/src/general/corrections.jl @@ -458,6 +458,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 diff --git a/src/schemes/fluid/entropically_damped_sph/rhs.jl b/src/schemes/fluid/entropically_damped_sph/rhs.jl index 2ea1e4dc4d..99adb7352a 100644 --- a/src/schemes/fluid/entropically_damped_sph/rhs.jl +++ b/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, diff --git a/src/schemes/fluid/entropically_damped_sph/system.jl b/src/schemes/fluid/entropically_damped_sph/system.jl index 1ea4b5b74e..077a0f051a 100644 --- a/src/schemes/fluid/entropically_damped_sph/system.jl +++ b/src/schemes/fluid/entropically_damped_sph/system.jl @@ -120,9 +120,10 @@ function EntropicallyDampedSPHSystem(initial_condition; smoothing_kernel, smooth throw(ArgumentError("`acceleration` must be of length $NDIMS for a $(NDIMS)D problem")) end + check_akinci_correction(surface_tension, correction) + 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 a surface-normal method")) end @@ -253,6 +254,12 @@ end system_correction(system::EntropicallyDampedSPHSystem) = system.correction +@inline function surface_normal_density(system::EntropicallyDampedSPHSystem, particle, + density) + return surface_normal_density(system, system.surface_tension, system.correction, + system.density_calculator, particle, density) +end + @inline function current_velocity(v, system::EntropicallyDampedSPHSystem) return view(v, 1:ndims(system), :) end @@ -299,11 +306,25 @@ function update_quantities!(system::EntropicallyDampedSPHSystem, v, u, end function update_pressure!(system::EntropicallyDampedSPHSystem, v, u, v_ode, u_ode, semi, t) + compute_akinci_correction_density!(system, system.correction, + system.density_calculator, u, u_ode, semi) + compute_surface_normal!(system, system.surface_normal_method, v, u, v_ode, u_ode, semi, t) compute_surface_delta_function!(system, system.surface_tension, semi) end +function compute_akinci_correction_density!(system, correction, density_calculator, + u, u_ode, semi) + return system +end + +function compute_akinci_correction_density!(system, ::AkinciFreeSurfaceCorrection, + density_calculator, u, u_ode, semi) + compute_akinci_correction_density!(system, density_calculator, u, u_ode, semi) + return system +end + function update_final!(system::EntropicallyDampedSPHSystem, v, u, v_ode, u_ode, semi, t; kwargs...) (; surface_tension) = system diff --git a/src/schemes/fluid/surface_normal_sph.jl b/src/schemes/fluid/surface_normal_sph.jl index 3562b93ddf..66f5207882 100644 --- a/src/schemes/fluid/surface_normal_sph.jl +++ b/src/schemes/fluid/surface_normal_sph.jl @@ -17,8 +17,13 @@ end function ColorfieldSurfaceNormal(; boundary_contact_threshold=0.1, interface_threshold=0.01, ideal_density_threshold=0.0) - thresholds = promote(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 @@ -47,6 +52,8 @@ end return extract_svector(cache.surface_normal, particle_system, particle) end +@inline surface_normal_density(system, particle, density) = density + function calc_normal!(system, neighbor_system, u_system, v, v_neighbor_system, u_neighbor_system, semi, surface_normal_method, neighbor_surface_normal_method) @@ -72,6 +79,8 @@ function calc_normal!(system::AbstractFluidSystem, neighbor_system::AbstractFlui m_b = hydrodynamic_mass(neighbor_system, neighbor) density_neighbor = current_density(v_neighbor_system, neighbor_system, neighbor) + density_neighbor = surface_normal_density(neighbor_system, neighbor, + density_neighbor) grad_kernel = smoothing_kernel_grad(system, pos_diff, distance, particle) for i in 1:ndims(system) cache.surface_normal[i, particle] += m_b / density_neighbor * grad_kernel[i] diff --git a/src/schemes/fluid/surface_tension.jl b/src/schemes/fluid/surface_tension.jl index 53e6151363..86f2a7b5dd 100644 --- a/src/schemes/fluid/surface_tension.jl +++ b/src/schemes/fluid/surface_tension.jl @@ -126,6 +126,15 @@ end @inline requires_surface_normal(::CohesionForceAkinci) = false @inline requires_surface_normal(::Any) = true +function check_akinci_correction(surface_tension, correction) + if surface_tension isa SurfaceTensionAkinci && + !(correction isa AkinciFreeSurfaceCorrection) + throw(ArgumentError("`SurfaceTensionAkinci` requires `AkinciFreeSurfaceCorrection`")) + end + + return surface_tension +end + function create_cache_surface_tension(::SurfaceTensionMomentumMorris, ELTYPE, NDIMS, nparticles) delta_s = Array{ELTYPE, 1}(undef, nparticles) @@ -145,21 +154,32 @@ end # By using the `@fastpow` macro, we are consciously trading off some precision in the result # for enhanced computational speed. This is especially useful in scenarios where performance # is a higher priority than exact precision. -@fastpow @inline function cohesion_force_akinci(surface_tension, support_radius, m_b, - pos_diff, distance) - (; surface_tension_coefficient) = surface_tension +@inline function cohesion_force_akinci(surface_tension, support_radius, m_b, pos_diff, + distance) + return cohesion_force_akinci(surface_tension.surface_tension_coefficient, + support_radius, m_b, pos_diff, distance) +end + +@fastpow @inline function cohesion_force_akinci(surface_tension_coefficient::Real, + support_radius, m_b, pos_diff, distance) + distance >= support_radius && return zero(pos_diff) - # Eq. 2 + # Eq. 2 in dimensionless form avoids scale-dependent powers up to `support_radius^9`. # We only reach this function when `sqrt(eps()) < distance <= support_radius` - if distance > 0.5 * support_radius + normalized_distance = distance / support_radius + if normalized_distance > one(normalized_distance) / 2 # Attractive force - C = (support_radius - distance)^3 * distance^3 + C = (1 - normalized_distance)^3 * normalized_distance^3 else - # `distance < 0.5 * support_radius` + # `distance <= 0.5 * support_radius` # Repulsive force - C = 2 * (support_radius - distance)^3 * distance^3 - support_radius^6 / 64.0 + C = 2 * (1 - normalized_distance)^3 * normalized_distance^3 - + one(normalized_distance) / 64 end - C *= 32.0 / (pi * support_radius^9) + normalization = oftype(support_radius, 32 / pi) + normalization = ((normalization / support_radius) / support_radius) / + support_radius + C *= normalization # Eq. 1 in acceleration form cohesion_force = -surface_tension_coefficient * m_b * C * pos_diff / distance @@ -167,18 +187,36 @@ end return cohesion_force end +@inline function pair_support_radius_akinci(particle_system, neighbor_system, particle, + neighbor) + support_radius_a = compact_support(system_smoothing_kernel(particle_system), + smoothing_length(particle_system, particle)) + support_radius_b = compact_support(system_smoothing_kernel(neighbor_system), + smoothing_length(neighbor_system, neighbor)) + # Both directed interaction passes are guaranteed to contain pairs inside this radius. + return min(support_radius_a, support_radius_b) +end + +# This symmetric numerical default preserves pair momentum. A physical multiphase model can +# specialize this function once interfacial coefficients are represented explicitly. +@inline function pair_surface_tension_coefficient_akinci(surface_tension_a, + surface_tension_b) + return (surface_tension_a.surface_tension_coefficient + + surface_tension_b.surface_tension_coefficient) / 2 +end + @inline function adhesion_force_akinci(surface_tension, support_radius, m_b, pos_diff, distance, adhesion_coefficient) distance >= support_radius && return zero(pos_diff) distance <= 0.5 * support_radius && return zero(pos_diff) - # Eq. 7. The factored radicand avoids cancellation close to the support boundary. - radicand = 2 * (2 * distance - support_radius) * - (support_radius - distance) / support_radius + # Eq. 7 in dimensionless form avoids cancellation and scale-dependent intermediates. + normalized_distance = distance / support_radius + radicand = 2 * (2 * normalized_distance - 1) * (1 - normalized_distance) fourth_root = sqrt(sqrt(max(zero(radicand), radicand))) - normalization = convert(typeof(support_radius), 0.007) / - (support_radius^3 * sqrt(sqrt(support_radius))) + normalization = convert(typeof(support_radius), 0.007) + normalization = ((normalization / support_radius) / support_radius) / support_radius A = normalization * fourth_root # Eq. 6 in acceleration form with `m_b` being the boundary mass calculated as @@ -204,17 +242,17 @@ end particle, neighbor, pos_diff, distance, rho_a, rho_b, grad_kernel, surface_tension_correction) - (; smoothing_kernel) = particle_system - # No cohesion with oneself. See `src/general/smoothing_kernels.jl` for more details. distance^2 < eps(initial_smoothing_length(particle_system)^2) && return dv_particle m_b = hydrodynamic_mass(neighbor_system, neighbor) - support_radius = compact_support(smoothing_kernel, - smoothing_length(particle_system, particle)) + support_radius = pair_support_radius_akinci(particle_system, neighbor_system, particle, + neighbor) + surface_tension_coefficient = pair_surface_tension_coefficient_akinci(surface_tension_a, + surface_tension_b) dv_particle[] += surface_tension_correction * - cohesion_force_akinci(surface_tension_a, support_radius, m_b, + cohesion_force_akinci(surface_tension_coefficient, support_radius, m_b, pos_diff, distance) return dv_particle @@ -228,23 +266,28 @@ end neighbor, pos_diff, distance, rho_a, rho_b, grad_kernel, surface_tension_correction) - (; smoothing_kernel) = particle_system - (; surface_tension_coefficient) = surface_tension_a - - smoothing_length_ = smoothing_length(particle_system, particle) # No surface tension with oneself. See `src/general/smoothing_kernels.jl` for more details. distance^2 < eps(initial_smoothing_length(particle_system)^2) && return dv_particle + m_a = hydrodynamic_mass(particle_system, particle) m_b = hydrodynamic_mass(neighbor_system, neighbor) n_a = surface_normal(particle_system, particle) n_b = surface_normal(neighbor_system, neighbor) - support_radius = compact_support(smoothing_kernel, smoothing_length_) + smoothing_length_a = smoothing_length(particle_system, particle) + smoothing_length_b = smoothing_length(neighbor_system, neighbor) + support_radius = pair_support_radius_akinci(particle_system, neighbor_system, particle, + neighbor) + distance >= support_radius && return dv_particle + surface_tension_coefficient = pair_surface_tension_coefficient_akinci(surface_tension_a, + surface_tension_b) dv_particle[] += surface_tension_correction * - cohesion_force_akinci(surface_tension_a, support_radius, m_b, + cohesion_force_akinci(surface_tension_coefficient, support_radius, m_b, pos_diff, distance) + normal_difference = smoothing_length_a * n_a - smoothing_length_b * n_b + neighbor_mass_weight = 2 * m_b / (m_a + m_b) dv_particle[] -= surface_tension_correction * surface_tension_coefficient * - (n_a - n_b) * smoothing_length_ + neighbor_mass_weight * normal_difference return dv_particle end diff --git a/src/schemes/fluid/weakly_compressible_sph/rhs.jl b/src/schemes/fluid/weakly_compressible_sph/rhs.jl index 8ceb29d27f..5cbc31fd60 100644 --- a/src/schemes/fluid/weakly_compressible_sph/rhs.jl +++ b/src/schemes/fluid/weakly_compressible_sph/rhs.jl @@ -71,10 +71,15 @@ function interact!(dv, v_particle_system, u_particle_system, # Determine correction factors. # This can usually be ignored, as these are all 1 when no correction is used. + 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, - rho_a, rho_b) + correction_rho_a, + correction_rho_b) # For `ContinuityDensity` without correction, this is equivalent to # dv_pressure = -m_b * (p_a + p_b) / (rho_a * rho_b) * grad_kernel @@ -126,6 +131,20 @@ function interact!(dv, v_particle_system, u_particle_system, return dv end +@inline function correction_density(::AkinciFreeSurfaceCorrection, + system::Union{WeaklyCompressibleSPHSystem, + EntropicallyDampedSPHSystem}, + particle, density) + if system.density_calculator isa ContinuityDensity && + haskey(system.cache, :kernel_summation_density) + return @inbounds system.cache.kernel_summation_density[particle] + end + + return density +end + +@inline correction_density(correction, system, particle, density) = density + @propagate_inbounds function neighbor_pressure(v_neighbor_system, neighbor_system, neighbor, p_a) return current_pressure(v_neighbor_system, neighbor_system, neighbor) diff --git a/src/schemes/fluid/weakly_compressible_sph/system.jl b/src/schemes/fluid/weakly_compressible_sph/system.jl index c97c829824..5172b50759 100644 --- a/src/schemes/fluid/weakly_compressible_sph/system.jl +++ b/src/schemes/fluid/weakly_compressible_sph/system.jl @@ -131,6 +131,8 @@ function WeaklyCompressibleSPHSystem(initial_condition; smoothing_kernel, throw(ArgumentError("`ShepardKernelCorrection` cannot be used with `ContinuityDensity`")) end + check_akinci_correction(surface_tension, correction) + surface_normal_method = default_surface_normal_method(surface_tension, surface_normal_method) @@ -245,6 +247,23 @@ end system_correction(system::WeaklyCompressibleSPHSystem) = system.correction +@inline function surface_normal_density(system::WeaklyCompressibleSPHSystem, particle, + density) + return surface_normal_density(system, system.surface_tension, system.correction, + system.density_calculator, particle, density) +end + +@inline function surface_normal_density(system, surface_tension, correction, + density_calculator, particle, density) + return density +end + +@inline function surface_normal_density(system, ::SurfaceTensionAkinci, + ::AkinciFreeSurfaceCorrection, + ::ContinuityDensity, particle, density) + return @inbounds system.cache.kernel_summation_density[particle] +end + @propagate_inbounds function current_velocity(v, system::WeaklyCompressibleSPHSystem) return current_velocity(v, system.density_calculator, system) end @@ -338,6 +357,22 @@ function update_pressure!(system::WeaklyCompressibleSPHSystem, v, u, v_ode, u_od return system end +function compute_correction_values!(system::WeaklyCompressibleSPHSystem, + ::AkinciFreeSurfaceCorrection, u, + v_ode, u_ode, semi) + compute_akinci_correction_density!(system, system.density_calculator, u, u_ode, semi) + return system +end + +function compute_akinci_correction_density!(system, ::ContinuityDensity, u, u_ode, semi) + summation_density!(system, semi, u, u_ode, system.cache.kernel_summation_density) + return system +end + +function compute_akinci_correction_density!(system, ::SummationDensity, u, u_ode, semi) + return system +end + function update_final!(system::WeaklyCompressibleSPHSystem, v, u, v_ode, u_ode, semi, t; kwargs...) (; surface_tension) = system diff --git a/test/schemes/fluid/surface_normal_sph.jl b/test/schemes/fluid/surface_normal_sph.jl index 5eb8a81704..a4be877616 100644 --- a/test/schemes/fluid/surface_normal_sph.jl +++ b/test/schemes/fluid/surface_normal_sph.jl @@ -76,6 +76,10 @@ function create_fluid_system(coordinates, velocity, mass, density, particle_spac system = WeaklyCompressibleSPHSystem(fluid; smoothing_kernel, smoothing_length, density_calculator=SummationDensity(), state_equation, + correction=surface_tension isa + SurfaceTensionAkinci ? + AkinciFreeSurfaceCorrection(first(density)) : + nothing, surface_normal_method, reference_particle_spacing=particle_spacing, surface_tension) diff --git a/test/schemes/fluid/surface_tension.jl b/test/schemes/fluid/surface_tension.jl index d470da0ba9..33b5901c2a 100644 --- a/test/schemes/fluid/surface_tension.jl +++ b/test/schemes/fluid/surface_tension.jl @@ -27,6 +27,20 @@ interface_threshold=0.01f0, ideal_density_threshold=0.0f0) isa ColorfieldSurfaceNormal{Float32} + + invalid_thresholds = ((NaN, 0.01, 0.0), + (-Inf, 0.01, 0.0), + (0.1, Inf, 0.0), + (0.1, 0.01, 1.0im), + ("invalid", 0.01, 0.0)) + for (boundary_threshold, interface_threshold, density_threshold) in + invalid_thresholds + + @test_throws ArgumentError ColorfieldSurfaceNormal(; + boundary_contact_threshold=boundary_threshold, + interface_threshold, + ideal_density_threshold=density_threshold) + end end @testset "cohesion-only systems do not require normals" begin @@ -55,6 +69,14 @@ @test !haskey(system.cache, :surface_normal) @test !haskey(system.cache, :neighbor_count) @test !haskey(system.cache, :reference_particle_spacing) + + semi = Semidiscretization(system) + ode = semidiscretize(semi, (0.0, 0.1)) + v_ode, u_ode = ode.u0.x + dv_ode = zero(v_ode) + @test_nowarn TrixiParticles.kick!(dv_ode, v_ode, u_ode, ode.p, 0.0) + @test all(isfinite, dv_ode) + @test any(!iszero, dv_ode) end @test_throws ArgumentError WeaklyCompressibleSPHSystem(initial_condition; @@ -64,13 +86,32 @@ state_equation=StateEquationCole(sound_speed=10.0, reference_density=1.0, exponent=1), - surface_tension=SurfaceTensionAkinci()) + surface_tension=SurfaceTensionAkinci(), + correction=AkinciFreeSurfaceCorrection(1.0)) + @test_throws ArgumentError EntropicallyDampedSPHSystem(initial_condition; + smoothing_kernel, + smoothing_length, + sound_speed=10.0, + density_calculator=SummationDensity(), + surface_tension=SurfaceTensionAkinci(), + correction=AkinciFreeSurfaceCorrection(1.0)) + + @test_throws ArgumentError WeaklyCompressibleSPHSystem(initial_condition; + smoothing_kernel, + smoothing_length, + density_calculator=SummationDensity(), + state_equation=StateEquationCole(sound_speed=10.0, + reference_density=1.0, + exponent=1), + surface_tension=SurfaceTensionAkinci(), + reference_particle_spacing=1.0) @test_throws ArgumentError EntropicallyDampedSPHSystem(initial_condition; smoothing_kernel, smoothing_length, sound_speed=10.0, density_calculator=SummationDensity(), - surface_tension=SurfaceTensionAkinci()) + surface_tension=SurfaceTensionAkinci(), + reference_particle_spacing=1.0) full_akinci = WeaklyCompressibleSPHSystem(initial_condition; smoothing_kernel, smoothing_length, @@ -79,11 +120,47 @@ reference_density=1.0, exponent=1), surface_tension=SurfaceTensionAkinci(), + correction=AkinciFreeSurfaceCorrection(1.0), reference_particle_spacing=1.0) @test full_akinci.surface_normal_method isa ColorfieldSurfaceNormal @test haskey(full_akinci.cache, :surface_normal) end + @testset "EDAC applies the Akinci free-surface correction" begin + function initial_acceleration(density_calculator, correction) + rho0 = 1.0 + initial_condition = InitialCondition(; coordinates=[0.0 1.0; 0.0 0.0], + velocity=zeros(2, 2), mass=ones(2), + density=ones(2), particle_spacing=1.0) + system = EntropicallyDampedSPHSystem(initial_condition; + smoothing_kernel=WendlandC2Kernel{2}(), + smoothing_length=1.0, sound_speed=10.0, + density_calculator, correction, + surface_tension=CohesionForceAkinci(surface_tension_coefficient=0.1)) + semi = Semidiscretization(system) + ode = semidiscretize(semi, (0.0, 0.1)) + v_ode, u_ode = ode.u0.x + dv_ode = zero(v_ode) + TrixiParticles.kick!(dv_ode, v_ode, u_ode, ode.p, 0.0) + dv = TrixiParticles.wrap_v(dv_ode, system, semi) + v = TrixiParticles.wrap_v(v_ode, system, semi) + density = TrixiParticles.current_density(v, system, 1) + correction_density = TrixiParticles.correction_density(correction, system, 1, + density) + return copy(dv[1:2, :]), rho0 / correction_density + end + + for density_calculator in (SummationDensity(), ContinuityDensity()) + corrected, + correction_factor = initial_acceleration(density_calculator, + AkinciFreeSurfaceCorrection(1.0)) + uncorrected, _ = initial_acceleration(density_calculator, nothing) + @test correction_factor > 1 + @test maximum(abs, uncorrected) > 0 + @test isapprox(corrected, correction_factor * uncorrected; rtol=2eps()) + end + end + @testset "zero Morris coefficient does not restrict the time step" begin function calculate_initial_dt(surface_tension) initial_condition = InitialCondition(; coordinates=[0.0 1.0; 0.0 0.0], @@ -155,6 +232,20 @@ pos_diff, test_distance) * test_distance @test isapprox(zero[1], 0.0, atol=6e-15) @test isapprox(zero[2], 0.0, atol=6e-15) + + surface_tension_f32 = CohesionForceAkinci(surface_tension_coefficient=1.0f0) + for support_radius in (1.0f12, 1.0f-12) + distance = 0.75f0 * support_radius + force = TrixiParticles.cohesion_force_akinci(surface_tension_f32, + support_radius, 1.0f0, + Float32[distance, 0], distance) + expected = Float32(-32 / pi * (1 - 0.75)^3 * 0.75^3 / + Float64(support_radius)^3) + @test eltype(force) == Float32 + @test all(isfinite, force) + @test isapprox(force[1], expected; rtol=4eps(Float32)) + @test iszero(force[2]) + end end @testset verbose=true "adhesion_force_akinci" begin @@ -212,6 +303,76 @@ @test eltype(near_support) == Float32 @test all(isfinite, near_support) @test 0 < norm(near_support) < eps(Float32) + + for support_radius in (1.0f12, 1.0f-13) + distance = 0.75f0 * support_radius + force = TrixiParticles.adhesion_force_akinci(surface_tension, + support_radius, 1.0f0, + Float32[distance, 0], distance, + 1.0f0) + expected = Float32(-0.007 / Float64(support_radius)^3 / sqrt(2)) + @test all(isfinite, force) + @test isapprox(force[1], expected; rtol=4eps(Float32)) + @test iszero(force[2]) + end + end + + @testset "adaptive Akinci pair force conserves momentum" begin + surface_tension_a = SurfaceTensionAkinci(surface_tension_coefficient=0.8) + surface_tension_b = SurfaceTensionAkinci(surface_tension_coefficient=0.2) + correction = AkinciFreeSurfaceCorrection(1.0) + smoothing_kernel = WendlandC2Kernel{2}() + state_equation = StateEquationCole(sound_speed=10.0, reference_density=1.0, + exponent=1) + + function system_at(coordinate, mass, density, smoothing_length, surface_tension) + initial_condition = InitialCondition(; + coordinates=reshape([coordinate, 0.0], + 2, 1), + velocity=zeros(2, 1), mass=[mass], + density=[density], particle_spacing=0.25) + return WeaklyCompressibleSPHSystem(initial_condition; smoothing_kernel, + smoothing_length, + density_calculator=SummationDensity(), + state_equation, surface_tension, correction, + reference_particle_spacing=0.25) + end + + mass_a, mass_b = 1.0, 4.0 + smoothing_length_a, smoothing_length_b = 0.5, 1.0 + system_a = system_at(0.0, mass_a, 1.0, smoothing_length_a, + surface_tension_a) + system_b = system_at(0.6, mass_b, 1.0, smoothing_length_b, + surface_tension_b) + system_a.cache.surface_normal[:, 1] .= (1 / smoothing_length_a, 0.0) + system_b.cache.surface_normal[:, 1] .= (0.0, 0.0) + + pos_diff = SVector(-0.6, 0.0) + acceleration_a = Ref(zero(pos_diff)) + acceleration_b = Ref(zero(pos_diff)) + TrixiParticles.surface_tension_force!(acceleration_a, surface_tension_a, + surface_tension_b, system_a, system_b, + 1, 1, pos_diff, norm(pos_diff), 1.0, 1.0, + zero(pos_diff), 1) + TrixiParticles.surface_tension_force!(acceleration_b, surface_tension_b, + surface_tension_a, system_b, system_a, + 1, 1, -pos_diff, norm(pos_diff), 1.0, 1.0, + zero(pos_diff), 1) + + @test mass_a * acceleration_a[] ≈ -mass_b * acceleration_b[] + @test TrixiParticles.pair_surface_tension_coefficient_akinci(surface_tension_a, + surface_tension_b) == + 0.5 + @test TrixiParticles.pair_support_radius_akinci(system_a, system_b, 1, 1) ≈ 1.0 + + outside_shared_support = SVector(-1.2, 0.0) + acceleration = Ref(zero(outside_shared_support)) + TrixiParticles.surface_tension_force!(acceleration, surface_tension_a, + surface_tension_b, system_a, system_b, + 1, 1, outside_shared_support, + norm(outside_shared_support), 1.0, 1.0, + zero(outside_shared_support), 1) + @test iszero(acceleration[]) end @testset "compute_stress_tensors! (MomentumMorris)" begin diff --git a/test/systems/rigid_system.jl b/test/systems/rigid_system.jl index d549363d9d..ea0e32ca7b 100644 --- a/test/systems/rigid_system.jl +++ b/test/systems/rigid_system.jl @@ -500,6 +500,7 @@ density_calculator=SummationDensity(), state_equation, surface_tension=SurfaceTensionAkinci(surface_tension_coefficient=0.05), + correction=AkinciFreeSurfaceCorrection(fluid_density), reference_particle_spacing=particle_spacing) boundary_coordinates = reshape([1.5, 0.0], 2, 1) From cc239f46a707adaf27b178717fec5a84bfd42854 Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Fri, 14 Aug 2026 16:34:39 +0200 Subject: [PATCH 05/13] Scope Akinci core hardening fixes --- NEWS.md | 7 +- docs/src/systems/fluid.md | 41 +----- examples/fluid/falling_water_spheres_2d.jl | 3 - src/general/corrections.jl | 4 - .../fluid/entropically_damped_sph/rhs.jl | 15 +-- .../fluid/entropically_damped_sph/system.jl | 23 +--- src/schemes/fluid/surface_normal_sph.jl | 4 - src/schemes/fluid/surface_tension.jl | 66 +++------- .../fluid/weakly_compressible_sph/rhs.jl | 21 +--- .../fluid/weakly_compressible_sph/system.jl | 35 ------ test/schemes/fluid/surface_normal_sph.jl | 4 - test/schemes/fluid/surface_tension.jl | 117 +----------------- test/systems/rigid_system.jl | 1 - 13 files changed, 31 insertions(+), 310 deletions(-) diff --git a/NEWS.md b/NEWS.md index 49bc89f880..80fc3dc823 100644 --- a/NEWS.md +++ b/NEWS.md @@ -8,10 +8,9 @@ used in the Julia ecosystem. Notable changes will be documented in this file for ### Important Bugfixes -- Hardened surface tension model configuration by validating coefficients, avoiding unnecessary - surface-normal allocation for `CohesionForceAkinci`, requiring the published free-surface - correction for the full Akinci model, stabilizing Akinci kernels across floating-point scales, - and conserving pair momentum across unequal particle masses and smoothing lengths. +- 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 diff --git a/docs/src/systems/fluid.md b/docs/src/systems/fluid.md index b518f2de85..1e6289fbea 100644 --- a/docs/src/systems/fluid.md +++ b/docs/src/systems/fluid.md @@ -240,15 +240,13 @@ where: #### Normalization of surface normals -Models that require only the normal direction normalize the calculated gradient: +The calculated normals are normalized to unit vectors: ```math \hat{n}_a = \frac{n_a}{\Vert n_a \Vert}. ``` -The Morris models use this unit normal. The Akinci surface-area term instead retains the raw -gradient magnitude and multiplies it by the particle smoothing length, as required by the -published formulation. +Normalization ensures that the magnitude of the normals does not bias the curvature calculations or the resulting surface tension forces. #### Handling noise and errors in normal calculation @@ -301,8 +299,6 @@ the fluid-fluid surface force. Wall adhesion is controlled independently by the 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. -The full Akinci model also requires `AkinciFreeSurfaceCorrection`, which supplies the -particle-neighborhood correction in Equations 4 and 5 of Akinci et al. ### [Akinci-based intra-particle force surface tension and wall adhesion model](@id akinci_ipf) @@ -337,39 +333,14 @@ C(r)=\frac{32}{\pi h_c^9} #### Surface area minimization force -The surface area minimization force models curvature reduction. Let ``n_i`` denote the raw color -gradient and ``N_i=h_i n_i`` its dimensionless Akinci normal. For equal particle masses, the -published pair force is +The surface area minimization force models the curvature reduction effects, aligning particle motion to reduce the interface's total area. +It acts based on the difference in surface normals: ```math -F_{a\leftarrow b}^{\text{curvature}} = -\sigma m_a(N_a-N_b), +F_{\text{curvature}} = -\sigma (n_a - n_b), ``` -For equal particle masses and smoothing lengths, this is the published Akinci discretization. -For adaptive resolution, TrixiParticles forms the dimensionless normals -``N_i=h_i n_i`` and uses the symmetric pair force - -```math -F_{a\leftarrow b}^{\text{curvature}} = --\sigma\frac{2m_am_b}{m_a+m_b}(N_a-N_b). -``` - -For a shared correction factor, the corresponding accelerations are equal and opposite after -multiplication by particle mass, and reduce to the published expression when ``m_a=m_b``. The -cohesion and curvature forces use the smaller of the two compact-support radii. This guarantees -that both directed neighborhood searches contain every active pair. Coefficients from two -different materials are averaged only as a conservative numerical default; this is not a physical -multiphase interfacial-tension law. A future multiphase model must represent the independent pair -coefficient and the pair reference density explicitly. - -Both Akinci contributions are multiplied by the required neighborhood-deficiency correction - -```math -K_{ab}=\frac{2\rho_0}{\rho_a+\rho_b}, -``` - -implemented by `AkinciFreeSurfaceCorrection`. With `ContinuityDensity`, an auxiliary summation -density is reconstructed for this correction so that missing neighbors remain observable. +where ``n_a`` and ``n_b`` are the surface normals of the interacting particles. #### Wall adhesion force diff --git a/examples/fluid/falling_water_spheres_2d.jl b/examples/fluid/falling_water_spheres_2d.jl index 208c1b71bc..1f0016ac6d 100644 --- a/examples/fluid/falling_water_spheres_2d.jl +++ b/examples/fluid/falling_water_spheres_2d.jl @@ -61,8 +61,6 @@ 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, @@ -70,7 +68,6 @@ sphere_surface_tension = EntropicallyDampedSPHSystem(sphere1; 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, diff --git a/src/general/corrections.jl b/src/general/corrections.jl index f6f8fcdf9d..d97eeaf341 100644 --- a/src/general/corrections.jl +++ b/src/general/corrections.jl @@ -458,10 +458,6 @@ 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 diff --git a/src/schemes/fluid/entropically_damped_sph/rhs.jl b/src/schemes/fluid/entropically_damped_sph/rhs.jl index 99adb7352a..2ea1e4dc4d 100644 --- a/src/schemes/fluid/entropically_damped_sph/rhs.jl +++ b/src/schemes/fluid/entropically_damped_sph/rhs.jl @@ -42,14 +42,6 @@ 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) @@ -73,12 +65,12 @@ function interact!(dv, v_particle_system, u_particle_system, rho_b, pos_diff, distance, grad_kernel, correction) - dv_particle = Ref(pressure_correction * dv_pressure) + dv_particle = Ref(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, viscosity_correction) + v_a, v_b, grad_kernel) # Extra terms in the momentum equation when using a shifting technique @inbounds dv_shifting!(dv_particle, shifting_technique(particle_system), @@ -91,8 +83,7 @@ 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, - surface_tension_correction) + rho_a, rho_b, grad_kernel, 1) @inbounds adhesion_force!(dv_particle, surface_tension_a, particle_system, neighbor_system, diff --git a/src/schemes/fluid/entropically_damped_sph/system.jl b/src/schemes/fluid/entropically_damped_sph/system.jl index 077a0f051a..1ea4b5b74e 100644 --- a/src/schemes/fluid/entropically_damped_sph/system.jl +++ b/src/schemes/fluid/entropically_damped_sph/system.jl @@ -120,10 +120,9 @@ function EntropicallyDampedSPHSystem(initial_condition; smoothing_kernel, smooth throw(ArgumentError("`acceleration` must be of length $NDIMS for a $(NDIMS)D problem")) end - check_akinci_correction(surface_tension, correction) - 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 a surface-normal method")) end @@ -254,12 +253,6 @@ end system_correction(system::EntropicallyDampedSPHSystem) = system.correction -@inline function surface_normal_density(system::EntropicallyDampedSPHSystem, particle, - density) - return surface_normal_density(system, system.surface_tension, system.correction, - system.density_calculator, particle, density) -end - @inline function current_velocity(v, system::EntropicallyDampedSPHSystem) return view(v, 1:ndims(system), :) end @@ -306,25 +299,11 @@ function update_quantities!(system::EntropicallyDampedSPHSystem, v, u, end function update_pressure!(system::EntropicallyDampedSPHSystem, v, u, v_ode, u_ode, semi, t) - compute_akinci_correction_density!(system, system.correction, - system.density_calculator, u, u_ode, semi) - compute_surface_normal!(system, system.surface_normal_method, v, u, v_ode, u_ode, semi, t) compute_surface_delta_function!(system, system.surface_tension, semi) end -function compute_akinci_correction_density!(system, correction, density_calculator, - u, u_ode, semi) - return system -end - -function compute_akinci_correction_density!(system, ::AkinciFreeSurfaceCorrection, - density_calculator, u, u_ode, semi) - compute_akinci_correction_density!(system, density_calculator, u, u_ode, semi) - return system -end - function update_final!(system::EntropicallyDampedSPHSystem, v, u, v_ode, u_ode, semi, t; kwargs...) (; surface_tension) = system diff --git a/src/schemes/fluid/surface_normal_sph.jl b/src/schemes/fluid/surface_normal_sph.jl index 66f5207882..40824530c2 100644 --- a/src/schemes/fluid/surface_normal_sph.jl +++ b/src/schemes/fluid/surface_normal_sph.jl @@ -52,8 +52,6 @@ end return extract_svector(cache.surface_normal, particle_system, particle) end -@inline surface_normal_density(system, particle, density) = density - function calc_normal!(system, neighbor_system, u_system, v, v_neighbor_system, u_neighbor_system, semi, surface_normal_method, neighbor_surface_normal_method) @@ -79,8 +77,6 @@ function calc_normal!(system::AbstractFluidSystem, neighbor_system::AbstractFlui m_b = hydrodynamic_mass(neighbor_system, neighbor) density_neighbor = current_density(v_neighbor_system, neighbor_system, neighbor) - density_neighbor = surface_normal_density(neighbor_system, neighbor, - density_neighbor) grad_kernel = smoothing_kernel_grad(system, pos_diff, distance, particle) for i in 1:ndims(system) cache.surface_normal[i, particle] += m_b / density_neighbor * grad_kernel[i] diff --git a/src/schemes/fluid/surface_tension.jl b/src/schemes/fluid/surface_tension.jl index 86f2a7b5dd..fbb210e0c7 100644 --- a/src/schemes/fluid/surface_tension.jl +++ b/src/schemes/fluid/surface_tension.jl @@ -126,15 +126,6 @@ end @inline requires_surface_normal(::CohesionForceAkinci) = false @inline requires_surface_normal(::Any) = true -function check_akinci_correction(surface_tension, correction) - if surface_tension isa SurfaceTensionAkinci && - !(correction isa AkinciFreeSurfaceCorrection) - throw(ArgumentError("`SurfaceTensionAkinci` requires `AkinciFreeSurfaceCorrection`")) - end - - return surface_tension -end - function create_cache_surface_tension(::SurfaceTensionMomentumMorris, ELTYPE, NDIMS, nparticles) delta_s = Array{ELTYPE, 1}(undef, nparticles) @@ -154,14 +145,10 @@ end # By using the `@fastpow` macro, we are consciously trading off some precision in the result # for enhanced computational speed. This is especially useful in scenarios where performance # is a higher priority than exact precision. -@inline function cohesion_force_akinci(surface_tension, support_radius, m_b, pos_diff, - distance) - return cohesion_force_akinci(surface_tension.surface_tension_coefficient, - support_radius, m_b, pos_diff, distance) -end +@fastpow @inline function cohesion_force_akinci(surface_tension, support_radius, m_b, + pos_diff, distance) + (; surface_tension_coefficient) = surface_tension -@fastpow @inline function cohesion_force_akinci(surface_tension_coefficient::Real, - support_radius, m_b, pos_diff, distance) distance >= support_radius && return zero(pos_diff) # Eq. 2 in dimensionless form avoids scale-dependent powers up to `support_radius^9`. @@ -187,24 +174,6 @@ end return cohesion_force end -@inline function pair_support_radius_akinci(particle_system, neighbor_system, particle, - neighbor) - support_radius_a = compact_support(system_smoothing_kernel(particle_system), - smoothing_length(particle_system, particle)) - support_radius_b = compact_support(system_smoothing_kernel(neighbor_system), - smoothing_length(neighbor_system, neighbor)) - # Both directed interaction passes are guaranteed to contain pairs inside this radius. - return min(support_radius_a, support_radius_b) -end - -# This symmetric numerical default preserves pair momentum. A physical multiphase model can -# specialize this function once interfacial coefficients are represented explicitly. -@inline function pair_surface_tension_coefficient_akinci(surface_tension_a, - surface_tension_b) - return (surface_tension_a.surface_tension_coefficient + - surface_tension_b.surface_tension_coefficient) / 2 -end - @inline function adhesion_force_akinci(surface_tension, support_radius, m_b, pos_diff, distance, adhesion_coefficient) distance >= support_radius && return zero(pos_diff) @@ -242,17 +211,17 @@ end particle, neighbor, pos_diff, distance, rho_a, rho_b, grad_kernel, surface_tension_correction) + (; smoothing_kernel) = particle_system + # No cohesion with oneself. See `src/general/smoothing_kernels.jl` for more details. distance^2 < eps(initial_smoothing_length(particle_system)^2) && return dv_particle m_b = hydrodynamic_mass(neighbor_system, neighbor) - support_radius = pair_support_radius_akinci(particle_system, neighbor_system, particle, - neighbor) - surface_tension_coefficient = pair_surface_tension_coefficient_akinci(surface_tension_a, - surface_tension_b) + support_radius = compact_support(smoothing_kernel, + smoothing_length(particle_system, particle)) dv_particle[] += surface_tension_correction * - cohesion_force_akinci(surface_tension_coefficient, support_radius, m_b, + cohesion_force_akinci(surface_tension_a, support_radius, m_b, pos_diff, distance) return dv_particle @@ -266,28 +235,23 @@ end neighbor, pos_diff, distance, rho_a, rho_b, grad_kernel, surface_tension_correction) + (; smoothing_kernel) = particle_system + (; surface_tension_coefficient) = surface_tension_a + + smoothing_length_ = smoothing_length(particle_system, particle) # No surface tension with oneself. See `src/general/smoothing_kernels.jl` for more details. distance^2 < eps(initial_smoothing_length(particle_system)^2) && return dv_particle - m_a = hydrodynamic_mass(particle_system, particle) m_b = hydrodynamic_mass(neighbor_system, neighbor) n_a = surface_normal(particle_system, particle) n_b = surface_normal(neighbor_system, neighbor) - smoothing_length_a = smoothing_length(particle_system, particle) - smoothing_length_b = smoothing_length(neighbor_system, neighbor) - support_radius = pair_support_radius_akinci(particle_system, neighbor_system, particle, - neighbor) - distance >= support_radius && return dv_particle - surface_tension_coefficient = pair_surface_tension_coefficient_akinci(surface_tension_a, - surface_tension_b) + support_radius = compact_support(smoothing_kernel, smoothing_length_) dv_particle[] += surface_tension_correction * - cohesion_force_akinci(surface_tension_coefficient, support_radius, m_b, + cohesion_force_akinci(surface_tension_a, support_radius, m_b, pos_diff, distance) - normal_difference = smoothing_length_a * n_a - smoothing_length_b * n_b - neighbor_mass_weight = 2 * m_b / (m_a + m_b) dv_particle[] -= surface_tension_correction * surface_tension_coefficient * - neighbor_mass_weight * normal_difference + (n_a - n_b) * smoothing_length_ return dv_particle end diff --git a/src/schemes/fluid/weakly_compressible_sph/rhs.jl b/src/schemes/fluid/weakly_compressible_sph/rhs.jl index 5cbc31fd60..8ceb29d27f 100644 --- a/src/schemes/fluid/weakly_compressible_sph/rhs.jl +++ b/src/schemes/fluid/weakly_compressible_sph/rhs.jl @@ -71,15 +71,10 @@ function interact!(dv, v_particle_system, u_particle_system, # Determine correction factors. # This can usually be ignored, as these are all 1 when no correction is used. - 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) + rho_a, rho_b) # For `ContinuityDensity` without correction, this is equivalent to # dv_pressure = -m_b * (p_a + p_b) / (rho_a * rho_b) * grad_kernel @@ -131,20 +126,6 @@ function interact!(dv, v_particle_system, u_particle_system, return dv end -@inline function correction_density(::AkinciFreeSurfaceCorrection, - system::Union{WeaklyCompressibleSPHSystem, - EntropicallyDampedSPHSystem}, - particle, density) - if system.density_calculator isa ContinuityDensity && - haskey(system.cache, :kernel_summation_density) - return @inbounds system.cache.kernel_summation_density[particle] - end - - return density -end - -@inline correction_density(correction, system, particle, density) = density - @propagate_inbounds function neighbor_pressure(v_neighbor_system, neighbor_system, neighbor, p_a) return current_pressure(v_neighbor_system, neighbor_system, neighbor) diff --git a/src/schemes/fluid/weakly_compressible_sph/system.jl b/src/schemes/fluid/weakly_compressible_sph/system.jl index 5172b50759..c97c829824 100644 --- a/src/schemes/fluid/weakly_compressible_sph/system.jl +++ b/src/schemes/fluid/weakly_compressible_sph/system.jl @@ -131,8 +131,6 @@ function WeaklyCompressibleSPHSystem(initial_condition; smoothing_kernel, throw(ArgumentError("`ShepardKernelCorrection` cannot be used with `ContinuityDensity`")) end - check_akinci_correction(surface_tension, correction) - surface_normal_method = default_surface_normal_method(surface_tension, surface_normal_method) @@ -247,23 +245,6 @@ end system_correction(system::WeaklyCompressibleSPHSystem) = system.correction -@inline function surface_normal_density(system::WeaklyCompressibleSPHSystem, particle, - density) - return surface_normal_density(system, system.surface_tension, system.correction, - system.density_calculator, particle, density) -end - -@inline function surface_normal_density(system, surface_tension, correction, - density_calculator, particle, density) - return density -end - -@inline function surface_normal_density(system, ::SurfaceTensionAkinci, - ::AkinciFreeSurfaceCorrection, - ::ContinuityDensity, particle, density) - return @inbounds system.cache.kernel_summation_density[particle] -end - @propagate_inbounds function current_velocity(v, system::WeaklyCompressibleSPHSystem) return current_velocity(v, system.density_calculator, system) end @@ -357,22 +338,6 @@ function update_pressure!(system::WeaklyCompressibleSPHSystem, v, u, v_ode, u_od return system end -function compute_correction_values!(system::WeaklyCompressibleSPHSystem, - ::AkinciFreeSurfaceCorrection, u, - v_ode, u_ode, semi) - compute_akinci_correction_density!(system, system.density_calculator, u, u_ode, semi) - return system -end - -function compute_akinci_correction_density!(system, ::ContinuityDensity, u, u_ode, semi) - summation_density!(system, semi, u, u_ode, system.cache.kernel_summation_density) - return system -end - -function compute_akinci_correction_density!(system, ::SummationDensity, u, u_ode, semi) - return system -end - function update_final!(system::WeaklyCompressibleSPHSystem, v, u, v_ode, u_ode, semi, t; kwargs...) (; surface_tension) = system diff --git a/test/schemes/fluid/surface_normal_sph.jl b/test/schemes/fluid/surface_normal_sph.jl index a4be877616..5eb8a81704 100644 --- a/test/schemes/fluid/surface_normal_sph.jl +++ b/test/schemes/fluid/surface_normal_sph.jl @@ -76,10 +76,6 @@ function create_fluid_system(coordinates, velocity, mass, density, particle_spac system = WeaklyCompressibleSPHSystem(fluid; smoothing_kernel, smoothing_length, density_calculator=SummationDensity(), state_equation, - correction=surface_tension isa - SurfaceTensionAkinci ? - AkinciFreeSurfaceCorrection(first(density)) : - nothing, surface_normal_method, reference_particle_spacing=particle_spacing, surface_tension) diff --git a/test/schemes/fluid/surface_tension.jl b/test/schemes/fluid/surface_tension.jl index 33b5901c2a..36b4233870 100644 --- a/test/schemes/fluid/surface_tension.jl +++ b/test/schemes/fluid/surface_tension.jl @@ -86,32 +86,13 @@ state_equation=StateEquationCole(sound_speed=10.0, reference_density=1.0, exponent=1), - surface_tension=SurfaceTensionAkinci(), - correction=AkinciFreeSurfaceCorrection(1.0)) + surface_tension=SurfaceTensionAkinci()) @test_throws ArgumentError EntropicallyDampedSPHSystem(initial_condition; smoothing_kernel, smoothing_length, sound_speed=10.0, density_calculator=SummationDensity(), - surface_tension=SurfaceTensionAkinci(), - correction=AkinciFreeSurfaceCorrection(1.0)) - - @test_throws ArgumentError WeaklyCompressibleSPHSystem(initial_condition; - smoothing_kernel, - smoothing_length, - density_calculator=SummationDensity(), - state_equation=StateEquationCole(sound_speed=10.0, - reference_density=1.0, - exponent=1), - surface_tension=SurfaceTensionAkinci(), - reference_particle_spacing=1.0) - @test_throws ArgumentError EntropicallyDampedSPHSystem(initial_condition; - smoothing_kernel, - smoothing_length, - sound_speed=10.0, - density_calculator=SummationDensity(), - surface_tension=SurfaceTensionAkinci(), - reference_particle_spacing=1.0) + surface_tension=SurfaceTensionAkinci()) full_akinci = WeaklyCompressibleSPHSystem(initial_condition; smoothing_kernel, smoothing_length, @@ -120,47 +101,11 @@ reference_density=1.0, exponent=1), surface_tension=SurfaceTensionAkinci(), - correction=AkinciFreeSurfaceCorrection(1.0), reference_particle_spacing=1.0) @test full_akinci.surface_normal_method isa ColorfieldSurfaceNormal @test haskey(full_akinci.cache, :surface_normal) end - @testset "EDAC applies the Akinci free-surface correction" begin - function initial_acceleration(density_calculator, correction) - rho0 = 1.0 - initial_condition = InitialCondition(; coordinates=[0.0 1.0; 0.0 0.0], - velocity=zeros(2, 2), mass=ones(2), - density=ones(2), particle_spacing=1.0) - system = EntropicallyDampedSPHSystem(initial_condition; - smoothing_kernel=WendlandC2Kernel{2}(), - smoothing_length=1.0, sound_speed=10.0, - density_calculator, correction, - surface_tension=CohesionForceAkinci(surface_tension_coefficient=0.1)) - semi = Semidiscretization(system) - ode = semidiscretize(semi, (0.0, 0.1)) - v_ode, u_ode = ode.u0.x - dv_ode = zero(v_ode) - TrixiParticles.kick!(dv_ode, v_ode, u_ode, ode.p, 0.0) - dv = TrixiParticles.wrap_v(dv_ode, system, semi) - v = TrixiParticles.wrap_v(v_ode, system, semi) - density = TrixiParticles.current_density(v, system, 1) - correction_density = TrixiParticles.correction_density(correction, system, 1, - density) - return copy(dv[1:2, :]), rho0 / correction_density - end - - for density_calculator in (SummationDensity(), ContinuityDensity()) - corrected, - correction_factor = initial_acceleration(density_calculator, - AkinciFreeSurfaceCorrection(1.0)) - uncorrected, _ = initial_acceleration(density_calculator, nothing) - @test correction_factor > 1 - @test maximum(abs, uncorrected) > 0 - @test isapprox(corrected, correction_factor * uncorrected; rtol=2eps()) - end - end - @testset "zero Morris coefficient does not restrict the time step" begin function calculate_initial_dt(surface_tension) initial_condition = InitialCondition(; coordinates=[0.0 1.0; 0.0 0.0], @@ -317,64 +262,6 @@ end end - @testset "adaptive Akinci pair force conserves momentum" begin - surface_tension_a = SurfaceTensionAkinci(surface_tension_coefficient=0.8) - surface_tension_b = SurfaceTensionAkinci(surface_tension_coefficient=0.2) - correction = AkinciFreeSurfaceCorrection(1.0) - smoothing_kernel = WendlandC2Kernel{2}() - state_equation = StateEquationCole(sound_speed=10.0, reference_density=1.0, - exponent=1) - - function system_at(coordinate, mass, density, smoothing_length, surface_tension) - initial_condition = InitialCondition(; - coordinates=reshape([coordinate, 0.0], - 2, 1), - velocity=zeros(2, 1), mass=[mass], - density=[density], particle_spacing=0.25) - return WeaklyCompressibleSPHSystem(initial_condition; smoothing_kernel, - smoothing_length, - density_calculator=SummationDensity(), - state_equation, surface_tension, correction, - reference_particle_spacing=0.25) - end - - mass_a, mass_b = 1.0, 4.0 - smoothing_length_a, smoothing_length_b = 0.5, 1.0 - system_a = system_at(0.0, mass_a, 1.0, smoothing_length_a, - surface_tension_a) - system_b = system_at(0.6, mass_b, 1.0, smoothing_length_b, - surface_tension_b) - system_a.cache.surface_normal[:, 1] .= (1 / smoothing_length_a, 0.0) - system_b.cache.surface_normal[:, 1] .= (0.0, 0.0) - - pos_diff = SVector(-0.6, 0.0) - acceleration_a = Ref(zero(pos_diff)) - acceleration_b = Ref(zero(pos_diff)) - TrixiParticles.surface_tension_force!(acceleration_a, surface_tension_a, - surface_tension_b, system_a, system_b, - 1, 1, pos_diff, norm(pos_diff), 1.0, 1.0, - zero(pos_diff), 1) - TrixiParticles.surface_tension_force!(acceleration_b, surface_tension_b, - surface_tension_a, system_b, system_a, - 1, 1, -pos_diff, norm(pos_diff), 1.0, 1.0, - zero(pos_diff), 1) - - @test mass_a * acceleration_a[] ≈ -mass_b * acceleration_b[] - @test TrixiParticles.pair_surface_tension_coefficient_akinci(surface_tension_a, - surface_tension_b) == - 0.5 - @test TrixiParticles.pair_support_radius_akinci(system_a, system_b, 1, 1) ≈ 1.0 - - outside_shared_support = SVector(-1.2, 0.0) - acceleration = Ref(zero(outside_shared_support)) - TrixiParticles.surface_tension_force!(acceleration, surface_tension_a, - surface_tension_b, system_a, system_b, - 1, 1, outside_shared_support, - norm(outside_shared_support), 1.0, 1.0, - zero(outside_shared_support), 1) - @test iszero(acceleration[]) - end - @testset "compute_stress_tensors! (MomentumMorris)" begin # 1. Define Minimal Initial Condition with 2 Particles in 2D coords = [0.0 1.0; diff --git a/test/systems/rigid_system.jl b/test/systems/rigid_system.jl index ea0e32ca7b..d549363d9d 100644 --- a/test/systems/rigid_system.jl +++ b/test/systems/rigid_system.jl @@ -500,7 +500,6 @@ density_calculator=SummationDensity(), state_equation, surface_tension=SurfaceTensionAkinci(surface_tension_coefficient=0.05), - correction=AkinciFreeSurfaceCorrection(fluid_density), reference_particle_spacing=particle_spacing) boundary_coordinates = reshape([1.5, 0.0], 2, 1) From 02d16890f5a901555b41aee86ccb3c9494bd5006 Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Fri, 14 Aug 2026 13:49:50 +0200 Subject: [PATCH 06/13] Clarify surface normal scope and visualization --- docs/src/systems/fluid.md | 123 ++++++++++++++++-- src/general/semidiscretization.jl | 5 +- .../fluid/entropically_damped_sph/system.jl | 10 +- src/schemes/fluid/surface_normal_sph.jl | 4 +- .../fluid/weakly_compressible_sph/system.jl | 13 +- src/schemes/structure/rigid_body/system.jl | 4 +- test/schemes/fluid/surface_tension.jl | 10 +- 7 files changed, 138 insertions(+), 31 deletions(-) diff --git a/docs/src/systems/fluid.md b/docs/src/systems/fluid.md index 1e6289fbea..68abf4ba88 100644 --- a/docs/src/systems/fluid.md +++ b/docs/src/systems/fluid.md @@ -215,17 +215,20 @@ Pages = [joinpath("general", "corrections.jl")] ### Overview of surface normal calculation in SPH -Surface normals provide the directionality of forces acting at the fluid interface. They are -used by the full Akinci model and both Morris models, but not by the cohesion-only Akinci model. -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] diff --git a/src/general/semidiscretization.jl b/src/general/semidiscretization.jl index a9ed769dc6..29da1dcd85 100644 --- a/src/general/semidiscretization.jl +++ b/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 diff --git a/src/schemes/fluid/entropically_damped_sph/system.jl b/src/schemes/fluid/entropically_damped_sph/system.jl index 1ea4b5b74e..5556b13546 100644 --- a/src/schemes/fluid/entropically_damped_sph/system.jl +++ b/src/schemes/fluid/entropically_damped_sph/system.jl @@ -50,9 +50,11 @@ 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 the surface tension model requires normals) +- `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 is needed when using a surface-normal method. - `color_value`: Integer label used for calculation of surface normals. @@ -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) diff --git a/src/schemes/fluid/surface_normal_sph.jl b/src/schemes/fluid/surface_normal_sph.jl index 40824530c2..c7e423491a 100644 --- a/src/schemes/fluid/surface_normal_sph.jl +++ b/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. diff --git a/src/schemes/fluid/weakly_compressible_sph/system.jl b/src/schemes/fluid/weakly_compressible_sph/system.jl index c97c829824..4b0a1ff414 100644 --- a/src/schemes/fluid/weakly_compressible_sph/system.jl +++ b/src/schemes/fluid/weakly_compressible_sph/system.jl @@ -53,9 +53,11 @@ See [Weakly Compressible SPH](@ref wcsph) for more details on the method. 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 the surface tension model requires normals) +- `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 is needed when using a surface-normal method. - `color_value`: Integer label used for calculation of surface normals. @@ -332,8 +334,9 @@ function update_pressure!(system::WeaklyCompressibleSPHSystem, v, u, v_ode, u_od kernel_correct_density!(system, v, u, v_ode, u_ode, semi, correction, density_calculator) - # These are only computed when using surface tension + # Interface normals are independent geometric data and are computed whenever a method is set. compute_surface_normal!(system, surface_normal_method, v, u, v_ode, u_ode, semi, t) + # The surface delta function is specific to surface-tension formulations that require it. compute_surface_delta_function!(system, surface_tension, semi) return system end @@ -342,7 +345,7 @@ function update_final!(system::WeaklyCompressibleSPHSystem, 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_shifting!(system, shifting_technique(system), v, u, v_ode, u_ode, semi) diff --git a/src/schemes/structure/rigid_body/system.jl b/src/schemes/structure/rigid_body/system.jl index dccacc2668..da2a93b609 100644 --- a/src/schemes/structure/rigid_body/system.jl +++ b/src/schemes/structure/rigid_body/system.jl @@ -36,7 +36,7 @@ torque and applied consistently to all rigid particles. colorfield initialization so fluids using [`ColorfieldSurfaceNormal`](@ref) can detect contact with rigid bodies, it participates in the multi-system color sanity check for - surface-tension setups, and it is written to VTK output as `"color"`. + colorfield-normal calculations, and it is written to VTK output as `"color"`. """ struct RigidBodySystem{BM, CTM, NDIMS, ELTYPE <: Real, IC, ARRAY1D, ARRAY2D, ST, CM, CMV, I, II, AV, RF, RT, AAF, GA, C} <: @@ -642,7 +642,7 @@ function check_configuration(system::RigidBodySystem, systems, nhs) throw(ArgumentError("`RigidBodySystem` with `BoundaryModelDummyParticles` " * "requires `reference_particle_spacing` to be set on " * "the boundary model when used together with " * - "`ColorfieldSurfaceNormal` or a surface tension model.")) + "`ColorfieldSurfaceNormal`.")) end end end diff --git a/test/schemes/fluid/surface_tension.jl b/test/schemes/fluid/surface_tension.jl index 36b4233870..687de166ef 100644 --- a/test/schemes/fluid/surface_tension.jl +++ b/test/schemes/fluid/surface_tension.jl @@ -58,11 +58,11 @@ state_equation=StateEquationCole(sound_speed=10.0, reference_density=1.0, exponent=1), - surface_tension) + surface_tension, color_value=0) edac = EntropicallyDampedSPHSystem(initial_condition; smoothing_kernel, smoothing_length, sound_speed=10.0, density_calculator=SummationDensity(), - surface_tension) + surface_tension, color_value=0) for system in (wcsph, edac) @test isnothing(system.surface_normal_method) @@ -79,6 +79,8 @@ @test any(!iszero, dv_ode) end + @test isnothing(TrixiParticles.check_system_color((wcsph, edac))) + @test_throws ArgumentError WeaklyCompressibleSPHSystem(initial_condition; smoothing_kernel, smoothing_length, @@ -101,9 +103,11 @@ reference_density=1.0, exponent=1), surface_tension=SurfaceTensionAkinci(), - reference_particle_spacing=1.0) + reference_particle_spacing=1.0, + color_value=0) @test full_akinci.surface_normal_method isa ColorfieldSurfaceNormal @test haskey(full_akinci.cache, :surface_normal) + @test_throws ArgumentError TrixiParticles.check_system_color((full_akinci, wcsph)) end @testset "zero Morris coefficient does not restrict the time step" begin From 300aebe7e96648320aeac6f959c4039375e1be88 Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Fri, 14 Aug 2026 14:15:59 +0200 Subject: [PATCH 07/13] Implement multicolor surface normals and threshold handling --- NEWS.md | 2 + docs/src/systems/fluid.md | 147 ++++++++++++++++-- src/general/semidiscretization.jl | 33 ++-- src/preprocessing/particle_packing/system.jl | 14 ++ src/schemes/boundary/wall_boundary/system.jl | 7 +- .../fluid/entropically_damped_sph/system.jl | 11 +- src/schemes/fluid/fluid.jl | 5 +- .../implicit_incompressible_sph/system.jl | 9 +- src/schemes/fluid/surface_normal_sph.jl | 130 +++++++++++----- src/schemes/fluid/surface_tension.jl | 11 +- .../fluid/weakly_compressible_sph/system.jl | 11 +- src/schemes/structure/rigid_body/system.jl | 4 +- test/schemes/fluid/surface_normal_sph.jl | 111 ++++++++++++- test/schemes/fluid/surface_tension.jl | 44 ++++-- test/systems/iisph_system.jl | 4 +- test/systems/packing_system.jl | 3 + 16 files changed, 429 insertions(+), 117 deletions(-) diff --git a/NEWS.md b/NEWS.md index 80fc3dc823..74dcaf1ffd 100644 --- a/NEWS.md +++ b/NEWS.md @@ -11,6 +11,8 @@ used in the Julia ecosystem. Notable changes will be documented in this file for - 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 diff --git a/docs/src/systems/fluid.md b/docs/src/systems/fluid.md index 68abf4ba88..3f0490af0f 100644 --- a/docs/src/systems/fluid.md +++ b/docs/src/systems/fluid.md @@ -241,6 +241,11 @@ 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``. +TrixiParticles evaluates this sum over every interacting physical fluid system. A neighboring +fluid therefore contributes its `color_value` even when it does not compute its own normals. +Particle-packing preprocessing systems are excluded. At a free surface, particles in the +unrepresented exterior phase are absent from the sum. + ```@eval using CairoMakie @@ -278,6 +283,118 @@ The color field is approximately constant within either phase. Its gradient is l 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. +##### Multiple color values + +With more than two color values, every transition between unequal values contributes to the +color-field gradient. The direction of each gradient points toward the larger color value, while +its magnitude depends on the size of the color jump. + +This is useful when several represented fluid phases or materials need to remain distinguishable in a +single scalar field. The value of ``c`` identifies the local region, while ``\nabla c`` locates and +orients each interface. This information can support interface reconstruction, phase-specific +boundary conditions, and post-processing. Assigning the same color value to two adjacent regions +deliberately makes their common boundary invisible to the color-field gradient. + +```@eval +using CairoMakie + +let + coordinate = range(-3.0, 3.0, length=601) + interface_width = 0.18 + smooth_step(position) = @. 0.5 * + (1.0 + tanh((coordinate - position) / interface_width)) + smooth_step_gradient(position) = @. 0.5 / interface_width / + cosh((coordinate - position) / interface_width)^2 + + color_a, color_b, color_c = 0.0, 2.0, 1.0 + colorfield = color_a .+ + (color_b - color_a) .* smooth_step(-1.0) .+ + (color_c - color_b) .* smooth_step(1.0) + colorfield_gradient = (color_b - color_a) .* smooth_step_gradient(-1.0) .+ + (color_c - color_b) .* smooth_step_gradient(1.0) + + fig = Figure(size=(1000, 430), fontsize=18) + color_axis = Axis(fig[1, 1], + xlabel="signed distance s/h", ylabel="color field c", + title="Three color values") + gradient_axis = Axis(fig[1, 2], + xlabel="signed distance s/h", ylabel="dc/d(s/h)", + title="Interface gradients") + + lines!(color_axis, coordinate, colorfield, color=:steelblue, linewidth=3) + lines!(gradient_axis, coordinate, colorfield_gradient, color=:darkorange, + linewidth=3) + text!(color_axis, -2.0, 0.3, text="A: c = 0", align=(:center, :center)) + text!(color_axis, 0.0, 1.5, text="B: c = 2", align=(:center, :center)) + text!(color_axis, 2.0, 0.7, text="C: c = 1", align=(:center, :center)) + vlines!(color_axis, [-1.0, 1.0], color=:black, linestyle=:dash, linewidth=2) + vlines!(gradient_axis, [-1.0, 1.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("multiple_color_values.png", fig) + + vertical_coordinate = range(-1.0, 1.0, length=101) + colorfield_2d = repeat(reshape(colorfield, :, 1), 1, length(vertical_coordinate)) + particle_coordinates = [(x, y) for x in -2.75:0.25:2.75 + for y in -0.8:0.25:0.8] + particle_x = first.(particle_coordinates) + particle_y = last.(particle_coordinates) + + normal_y = collect(range(-0.7, 0.7, length=5)) + normal_origins_x = vcat(fill(-1.0, length(normal_y)), + fill(1.0, length(normal_y))) + normal_origins_y = vcat(normal_y, normal_y) + normal_directions_x = vcat(fill(1.0, length(normal_y)), + fill(-1.0, length(normal_y))) + normal_directions_y = zeros(length(normal_directions_x)) + + normal_fig = Figure(size=(1000, 430), fontsize=18) + normal_axis = Axis(normal_fig[1, 1], aspect=DataAspect(), + xlabel="x/h", ylabel="y/h", + title="Unit normals at multiple color interfaces") + heatmap = heatmap!(normal_axis, coordinate, vertical_coordinate, colorfield_2d, + colormap=:viridis, colorrange=(0.0, 2.0)) + scatter!(normal_axis, particle_x, particle_y, color=(:black, 0.3), markersize=5) + text!(normal_axis, -2.0, 0.88, text="A: c = 0", color=:white, + align=(:center, :top)) + text!(normal_axis, 0.0, 0.88, text="B: c = 2", color=:black, + align=(:center, :top)) + text!(normal_axis, 2.0, 0.88, text="C: c = 1", color=:white, + align=(:center, :top)) + vlines!(normal_axis, [-1.0, 1.0], color=:white, linewidth=3) + arrows2d!(normal_axis, normal_origins_x, normal_origins_y, + normal_directions_x, normal_directions_y, + normalize=true, lengthscale=0.45, color=:black, + shaftwidth=3, tipwidth=14, tiplength=10) + xlims!(normal_axis, extrema(coordinate)) + ylims!(normal_axis, extrema(vertical_coordinate)) + Colorbar(normal_fig[1, 2], heatmap, label="color field c") + + CairoMakie.save("multiple_color_surface_normals.png", normal_fig) +end +``` + +![A color field and its gradient for three different color values](multiple_color_values.png) + +Here the regions from left to right have ``c_A=0``, ``c_B=2``, and ``c_C=1``. At the +``A-B`` interface, the positive gradient points from ``A`` to ``B``. At the ``B-C`` interface, +the negative gradient points from ``C`` to ``B``. The first peak is twice as large because +``|c_B-c_A|=2`` instead of ``|c_C-c_B|=1``. Thus, numerical color differences affect an +unnormalized color-field normal. Normalization removes this difference from the direction but not +from formulations that retain the raw gradient magnitude, including the Akinci area term and the +momentum-conserving Morris surface delta. Color contrasts must therefore be selected consistently +when those models are used. + +![Resulting surface normals for three color values](multiple_color_surface_normals.png) + +The unit-normal directions are defined only in the transition regions. At the ``A-B`` interface +they point to the right, from ``A`` toward the larger value in ``B``. At the ``B-C`` interface they +point to the left, again toward ``B``. Inside each constant-color region the gradient vanishes, so +the unit normal is undefined and is represented numerically by a zero vector. If the opposite +orientation is required, the normal sign must be reversed according to the chosen phase convention. + #### Normalization of surface normals The color-field gradient ``n_a`` is generally not a unit vector. Formulations that require only @@ -288,8 +405,10 @@ the interface orientation use the unit normal ``` 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. +gradient. In TrixiParticles, standalone analysis/VTK output and the Akinci surface-area force use +the filtered, unnormalized gradient. The Morris formulations use unit normals for curvature or +surface-stress calculations; the momentum-conserving formulation separately retains the raw +gradient magnitude as its surface delta function. ```@eval using CairoMakie @@ -346,12 +465,18 @@ without changing the interface geometry. #### Handling noise and errors in normal calculation 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: +incomplete kernel support can produce small or poorly resolved normal estimates. The +[`ColorfieldSurfaceNormal`](@ref) thresholds mitigate this as follows: 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. +2. `interface_threshold` rejects gradients for which the dimensionless magnitude + ``R\lVert n\rVert`` does not exceed the configured cutoff, where ``R`` is the kernel support + radius. This applies to standalone, Akinci, and Morris normals. +3. `ideal_density_threshold` optionally suppresses particles whose neighbor count is close to + ideal full support. This heuristic is intended only for a free surface with an unrepresented + exterior phase. It must remain zero for fully represented multiphase interfaces, where valid + interface particles can have full support. +4. Curvature calculations use a corrected formulation to reduce errors near interface fringes. ```@autodocs Modules = [TrixiParticles] @@ -430,14 +555,16 @@ C(r)=\frac{32}{\pi h_c^9} #### Surface area minimization force -The surface area minimization force models the curvature reduction effects, aligning particle motion to reduce the interface's total area. -It acts based on the difference in surface normals: +The surface area minimization term models curvature reduction by using the difference between the +raw color gradients. In the implementation it is evaluated in acceleration form as ```math -F_{\text{curvature}} = -\sigma (n_a - n_b), +a_{a,\text{area}} = -\sigma h_a (n_a - n_b), ``` -where ``n_a`` and ``n_b`` are the surface normals of the interacting particles. +where ``n_a`` and ``n_b`` are the unnormalized color gradients of the interacting particles and +``h_a`` is the smoothing length of particle ``a``. The factor ``h_a`` makes the color-normal term +dimensionless, consistent with the Akinci formulation. #### Wall adhesion force diff --git a/src/general/semidiscretization.jl b/src/general/semidiscretization.jl index 29da1dcd85..a0572ce75c 100644 --- a/src/general/semidiscretization.jl +++ b/src/general/semidiscretization.jl @@ -917,25 +917,30 @@ 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 system.surface_normal_method isa ColorfieldSurfaceNormal + return surface_normal_method(system) isa ColorfieldSurfaceNormal end if requires_color_check - - # 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 + 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 + + 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 diff --git a/src/preprocessing/particle_packing/system.jl b/src/preprocessing/particle_packing/system.jl index a4bb7e9aa7..3cb625677c 100644 --- a/src/preprocessing/particle_packing/system.jl +++ b/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)}(), diff --git a/src/schemes/boundary/wall_boundary/system.jl b/src/schemes/boundary/wall_boundary/system.jl index f9864ecc05..c6111e821d 100644 --- a/src/schemes/boundary/wall_boundary/system.jl +++ b/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} diff --git a/src/schemes/fluid/entropically_damped_sph/system.jl b/src/schemes/fluid/entropically_damped_sph/system.jl index 5556b13546..028b08090f 100644 --- a/src/schemes/fluid/entropically_damped_sph/system.jl +++ b/src/schemes/fluid/entropically_damped_sph/system.jl @@ -55,12 +55,11 @@ See [Entropically Damped Artificial Compressibility for SPH](@ref edac) for more 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 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 - reads the resulting boundary colorfield to detect wall contact. +- `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, diff --git a/src/schemes/fluid/fluid.jl b/src/schemes/fluid/fluid.jl index d94943f7bc..83bdae16ac 100644 --- a/src/schemes/fluid/fluid.jl +++ b/src/schemes/fluid/fluid.jl @@ -257,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), @@ -300,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) && diff --git a/src/schemes/fluid/implicit_incompressible_sph/system.jl b/src/schemes/fluid/implicit_incompressible_sph/system.jl index e4ef03584f..db5f5d7eb0 100644 --- a/src/schemes/fluid/implicit_incompressible_sph/system.jl +++ b/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, diff --git a/src/schemes/fluid/surface_normal_sph.jl b/src/schemes/fluid/surface_normal_sph.jl index c7e423491a..ceda09aa12 100644 --- a/src/schemes/fluid/surface_normal_sph.jl +++ b/src/schemes/fluid/surface_normal_sph.jl @@ -4,12 +4,24 @@ 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. +interface orientation. Every interacting fluid system contributes its `color_value` to the +discrete color gradient, even when that neighboring system does not compute its own normals. + +Without a surface-tension model and with [`SurfaceTensionAkinci`](@ref), the stored quantity is +the filtered, unnormalized color gradient. The Morris models store unit normals and retain the +raw gradient magnitude separately where required by the formulation. # Keywords -- `boundary_contact_threshold=0.1`: If this threshold is reached the fluid is assumed to be in contact with the boundary. -- `interface_threshold=0.01`: Threshold for normals to be removed as being invalid. -- `ideal_density_threshold=0.0`: Assume particles are inside if they are above this threshold, which is relative to the `ideal_neighbor_count`. +- `boundary_contact_threshold=0.1`: Finite value in `[0, 1]`. A dummy-boundary + particle is treated as being in contact with fluid when the magnitude of its smoothed + color field, normalized by the maximum magnitude, exceeds this value. +- `interface_threshold=0.01`: Finite, non-negative dimensionless cutoff ``\epsilon``. + A raw color gradient ``n`` is discarded when ``R\lVert n\rVert \leq \epsilon``, where + ``R`` is the kernel support radius. +- `ideal_density_threshold=0.0`: Finite value in `[0, 1]` controlling an optional + neighbor-count heuristic for free surfaces without a represented exterior phase. Zero + disables the heuristic. Keep this at zero for interfaces between represented phases, + since those interfaces can have full particle support. """ struct ColorfieldSurfaceNormal{ELTYPE} boundary_contact_threshold::ELTYPE @@ -19,14 +31,26 @@ end function ColorfieldSurfaceNormal(; boundary_contact_threshold=0.1, interface_threshold=0.01, ideal_density_threshold=0.0) - 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")) + boundary_threshold = validate_surface_normal_threshold(boundary_contact_threshold, + "boundary_contact_threshold"; + upper_bound=1) + normal_threshold = validate_surface_normal_threshold(interface_threshold, + "interface_threshold") + density_threshold = validate_surface_normal_threshold(ideal_density_threshold, + "ideal_density_threshold"; + upper_bound=1) + thresholds = promote(boundary_threshold, normal_threshold, density_threshold) + return ColorfieldSurfaceNormal(thresholds...) +end + +function validate_surface_normal_threshold(threshold, name; upper_bound=nothing) + if !(threshold isa Real) || !isfinite(threshold) || threshold < 0 || + (!isnothing(upper_bound) && threshold > upper_bound) + interval = isnothing(upper_bound) ? "non-negative" : "in [0, $upper_bound]" + throw(ArgumentError("`$name` must be a finite real number $interval")) end - thresholds = promote(thresholds...) - return ColorfieldSurfaceNormal(thresholds...) + return threshold end @inline function default_surface_normal_method(surface_tension, surface_normal_method) @@ -65,9 +89,11 @@ end # and Section 5 in Morris 2000 "Simulating surface tension with smoothed particle hydrodynamics". function calc_normal!(system::AbstractFluidSystem, neighbor_system::AbstractFluidSystem, u_system, v, - v_neighbor_system, u_neighbor_system, semi, surface_normal_method, - ::ColorfieldSurfaceNormal) + v_neighbor_system, u_neighbor_system, semi, + surface_normal_method::ColorfieldSurfaceNormal, + neighbor_surface_normal_method) (; cache) = system + color_b = neighbor_system.cache.color system_coords = current_coordinates(u_system, system) neighbor_system_coords = current_coordinates(u_neighbor_system, neighbor_system) @@ -81,7 +107,9 @@ function calc_normal!(system::AbstractFluidSystem, neighbor_system::AbstractFlui neighbor_system, neighbor) grad_kernel = smoothing_kernel_grad(system, pos_diff, distance, particle) for i in 1:ndims(system) - cache.surface_normal[i, particle] += m_b / density_neighbor * grad_kernel[i] + cache.surface_normal[i, + particle] += m_b / density_neighbor * color_b * + grad_kernel[i] end cache.neighbor_count[particle] += 1 @@ -106,19 +134,21 @@ function calc_boundary_normal!(system::AbstractFluidSystem, neighbor_system, u_s # TODO: move colorfield to extra step # TODO: this is only correct for a single fluid - # Reset to the constant boundary interpolated color values - colorfield .= initial_colorfield + # Contact detection depends on color magnitude, not interface orientation. + colorfield .= abs.(initial_colorfield) # Accumulate fluid neighbors foreach_point_neighbor(neighbor_system, system, neighbor_system_coords, system_coords, semi) do particle, neighbor, pos_diff, distance colorfield[particle] += hydrodynamic_mass(system, neighbor) / - current_density(v, system, neighbor) * system.cache.color * + current_density(v, system, neighbor) * + abs(system.cache.color) * smoothing_kernel(system, distance, particle) end maximum_colorfield = maximum(colorfield) + iszero(maximum_colorfield) && return system foreach_point_neighbor(system, neighbor_system, system_coords, neighbor_system_coords, @@ -148,43 +178,34 @@ end function remove_invalid_normals!(system::AbstractFluidSystem, surface_tension, surface_normal_method) - (; cache) = system - - # We remove invalid normals (too few neighbors) to reduce the impact of underdefined normals - for particle in each_integrated_particle(system) - # A corner has that many neighbors assuming a regular 2 * r distribution and a compact_support of 4r - if cache.neighbor_count[particle] < 2^ndims(system) + 1 - cache.surface_normal[1:ndims(system), particle] .= 0 - end - end - return system end -# See Morris 2000 "Simulating surface tension with smoothed particle hydrodynamics" function remove_invalid_normals!(system::AbstractFluidSystem, - surface_tension::Union{SurfaceTensionMorris, - SurfaceTensionMomentumMorris}, + surface_tension, surface_normal_method::ColorfieldSurfaceNormal) (; cache, smoothing_kernel) = system (; ideal_density_threshold, interface_threshold) = surface_normal_method (; neighbor_count) = cache smoothing_length_ = initial_smoothing_length(system) + support_radius = compact_support(smoothing_kernel, smoothing_length_) + minimum_neighbor_count = 2^ndims(system) + 1 - # We remove invalid normals i.e. they have a small norm (eq. 20) - normal_condition2 = (interface_threshold / - compact_support(smoothing_kernel, smoothing_length_))^2 + # Eq. 20 in Morris (2000) compares the color-gradient magnitude with ε/h. + normal_condition2 = (interface_threshold / support_radius)^2 + reset_surface_delta!(system, surface_tension) for particle in each_integrated_particle(system) - # Heuristic condition if there is no gas phase to find the free surface. - # We remove normals for particles which have a lot of support e.g. they are in the interior. - if ideal_density_threshold > 0 && - ideal_density_threshold * - ideal_neighbor_count(Val(ndims(system)), cache.reference_particle_spacing, - compact_support(smoothing_kernel, smoothing_length_)) < - neighbor_count[particle] + # This must stay disabled for fully supported interfaces between represented phases. + is_interior = ideal_density_threshold > 0 && + ideal_density_threshold * + ideal_neighbor_count(Val(ndims(system)), + cache.reference_particle_spacing, + support_radius) < neighbor_count[particle] + + if neighbor_count[particle] < minimum_neighbor_count || is_interior cache.surface_normal[1:ndims(system), particle] .= 0 continue end @@ -192,10 +213,15 @@ function remove_invalid_normals!(system::AbstractFluidSystem, particle_surface_normal = surface_normal(system, particle) norm2 = dot(particle_surface_normal, particle_surface_normal) - # See eq. 21 + # Eq. 21 in Morris (2000) defines the unit normal after rejecting weak gradients. if norm2 > normal_condition2 - cache.surface_normal[1:ndims(system), - particle] = particle_surface_normal / sqrt(norm2) + normal_magnitude = sqrt(norm2) + store_surface_delta!(system, surface_tension, particle, normal_magnitude) + + if normalize_surface_normals(surface_tension) + cache.surface_normal[1:ndims(system), + particle] = particle_surface_normal / normal_magnitude + end else cache.surface_normal[1:ndims(system), particle] .= 0 end @@ -204,6 +230,25 @@ function remove_invalid_normals!(system::AbstractFluidSystem, return system end +@inline normalize_surface_normals(surface_tension) = false +@inline normalize_surface_normals(::SurfaceTensionMorris) = true +@inline normalize_surface_normals(::SurfaceTensionMomentumMorris) = true + +@inline reset_surface_delta!(system, surface_tension) = system + +@inline function reset_surface_delta!(system, ::SurfaceTensionMomentumMorris) + set_zero!(system.cache.delta_s) + return system +end + +@inline store_surface_delta!(system, surface_tension, particle, normal_magnitude) = system + +@inline function store_surface_delta!(system, ::SurfaceTensionMomentumMorris, particle, + normal_magnitude) + system.cache.delta_s[particle] = normal_magnitude + return system +end + function compute_surface_normal!(system, surface_normal_method, v, u, v_ode, u_ode, semi, t) return system end @@ -217,7 +262,6 @@ function compute_surface_normal!(system::AbstractFluidSystem, set_zero!(cache.surface_normal) set_zero!(cache.neighbor_count) - # TODO: if color values are set only different systems need to be called @trixi_timeit timer() "compute surface normal" begin foreach_system_wrapped(semi, v_ode, u_ode) do neighbor_system, v_neighbor_system, diff --git a/src/schemes/fluid/surface_tension.jl b/src/schemes/fluid/surface_tension.jl index fbb210e0c7..8ea5b0ae95 100644 --- a/src/schemes/fluid/surface_tension.jl +++ b/src/schemes/fluid/surface_tension.jl @@ -317,16 +317,9 @@ function compute_surface_delta_function!(system, surface_tension, semi) return system end -# Eq. 6 in Morris 2000 "Simulating surface tension with smoothed particle hydrodynamics" +# The raw color-gradient magnitude is stored before normalization in +# `remove_invalid_normals!` (Eq. 6 in Morris 2000). function compute_surface_delta_function!(system, ::SurfaceTensionMomentumMorris, semi) - (; cache) = system - (; delta_s) = cache - - set_zero!(delta_s) - - @threaded semi for particle in each_integrated_particle(system) - delta_s[particle] = norm(surface_normal(system, particle)) - end return system end diff --git a/src/schemes/fluid/weakly_compressible_sph/system.jl b/src/schemes/fluid/weakly_compressible_sph/system.jl index 4b0a1ff414..76cdb9893e 100644 --- a/src/schemes/fluid/weakly_compressible_sph/system.jl +++ b/src/schemes/fluid/weakly_compressible_sph/system.jl @@ -58,12 +58,11 @@ See [Weakly Compressible SPH](@ref wcsph) for more details on the method. 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 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 - reads the resulting boundary colorfield to detect wall contact. +- `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 WeaklyCompressibleSPHSystem{NDIMS, ELTYPE <: Real, IC, MA, P, DC, SE, K, V, DD, COR, PF, SC, ST, B, SRFT, SRFN, PR, diff --git a/src/schemes/structure/rigid_body/system.jl b/src/schemes/structure/rigid_body/system.jl index da2a93b609..4f3b79090e 100644 --- a/src/schemes/structure/rigid_body/system.jl +++ b/src/schemes/structure/rigid_body/system.jl @@ -31,8 +31,8 @@ torque and applied consistently to all rigid particles. models when fluids interact with this rigid body. This is only evaluated for fluid-structure interaction with surface-tension-enabled fluid systems. -- `color_value`: Integer label stored as `system.cache.color`. - Currently this is used with `BoundaryModelDummyParticles` during +- `color_value`: Integer scalar stored as `system.cache.color`. + Its magnitude is used with `BoundaryModelDummyParticles` during colorfield initialization so fluids using [`ColorfieldSurfaceNormal`](@ref) can detect contact with rigid bodies, it participates in the multi-system color sanity check for diff --git a/test/schemes/fluid/surface_normal_sph.jl b/test/schemes/fluid/surface_normal_sph.jl index 5eb8a81704..099387eacb 100644 --- a/test/schemes/fluid/surface_normal_sph.jl +++ b/test/schemes/fluid/surface_normal_sph.jl @@ -62,6 +62,7 @@ end function create_fluid_system(coordinates, velocity, mass, density, particle_spacing, surface_tension; surface_normal_method=ColorfieldSurfaceNormal(), + color_value=1, NDIMS=2, smoothing_length=1.0, wall=false, walldistance=0.0, boundary_system_type=:wall, smoothing_kernel=SchoenbergCubicSplineKernel{NDIMS}()) @@ -78,7 +79,7 @@ function create_fluid_system(coordinates, velocity, mass, density, particle_spac state_equation, surface_normal_method, reference_particle_spacing=particle_spacing, - surface_tension) + surface_tension, color_value) if wall boundary_system = if boundary_system_type == :wall @@ -104,6 +105,114 @@ function create_fluid_system(coordinates, velocity, mass, density, particle_spac return system, boundary_system, semi, ode end +@testset "Multicolor fluid interfaces" begin + particle_spacing = 0.1 + smoothing_length = 0.15 + y_coordinates = collect(-0.5:particle_spacing:0.5) + coordinates_a = hcat(([x, y] for x in -0.5:particle_spacing:-0.1 + for y in y_coordinates)...) + coordinates_b = hcat(([x, y] for x in 0.0:particle_spacing:0.5 + for y in y_coordinates)...) + smoothing_kernel = WendlandC2Kernel{2}() + state_equation = StateEquationCole(sound_speed=10.0, reference_density=1000.0, + exponent=1) + normal_method = ColorfieldSurfaceNormal(interface_threshold=1.0e-6) + + function interface_normal(color_a, color_b) + initial_condition_a = InitialCondition(; coordinates=coordinates_a, + density=fill(1000.0, + size(coordinates_a, 2)), + particle_spacing) + initial_condition_b = InitialCondition(; coordinates=coordinates_b, + density=fill(1000.0, + size(coordinates_b, 2)), + particle_spacing) + + system_a = WeaklyCompressibleSPHSystem(initial_condition_a; smoothing_kernel, + smoothing_length, + density_calculator=SummationDensity(), + state_equation, + surface_normal_method=normal_method, + reference_particle_spacing=particle_spacing, + color_value=color_a) + # A fluid contributes its color even when it does not compute its own normals. + system_b = WeaklyCompressibleSPHSystem(initial_condition_b; smoothing_kernel, + smoothing_length, + density_calculator=SummationDensity(), + state_equation, color_value=color_b) + semi = Semidiscretization(system_a, system_b) + ode = semidiscretize(semi, (0.0, 0.01)) + v_ode, u_ode = ode.u0.x + TrixiParticles.update_systems_and_nhs(v_ode, u_ode, semi, 0.0) + + v_a = TrixiParticles.wrap_v(v_ode, system_a, semi) + u_a = TrixiParticles.wrap_u(u_ode, system_a, semi) + TrixiParticles.compute_surface_normal!(system_a, normal_method, v_a, u_a, + v_ode, u_ode, semi, 0.0) + + interface_particle = argmin(eachindex(eachcol(coordinates_a))) do particle + abs(coordinates_a[1, particle] + particle_spacing) + + abs(coordinates_a[2, particle]) + end + + return TrixiParticles.surface_normal(system_a, interface_particle), system_b + end + + increasing_normal, non_normal_neighbor = interface_normal(0, 2) + unit_jump_normal, _ = interface_normal(0, 1) + decreasing_normal, _ = interface_normal(2, 0) + equal_color_normal, _ = interface_normal(1, 1) + canceling_labels_normal, _ = interface_normal(-1, 1) + + @test isnothing(non_normal_neighbor.surface_normal_method) + @test increasing_normal[1] > 0 + @test decreasing_normal[1] < 0 + @test abs(increasing_normal[2]) < 100eps() + @test abs(decreasing_normal[2]) < 100eps() + @test isapprox(norm(increasing_normal), norm(decreasing_normal); rtol=1.0e-12) + @test isapprox(norm(increasing_normal), 2 * norm(unit_jump_normal); rtol=1.0e-12) + @test iszero(equal_color_normal) + @test canceling_labels_normal[1] > 0 +end + +@testset "Standalone surface-normal thresholds" begin + particle_spacing = 0.1 + coordinates = hcat(([x, y] for x in 0.0:particle_spacing:0.6 + for y in 0.0:particle_spacing:0.6)...) + velocity = zeros(2, size(coordinates, 2)) + density = fill(1000.0, size(coordinates, 2)) + mass = fill(10.0, size(coordinates, 2)) + + system, boundary, semi, + ode = create_fluid_system(coordinates, velocity, mass, density, particle_spacing, + nothing; + smoothing_length=0.15, + surface_normal_method=ColorfieldSurfaceNormal(interface_threshold=1.0e6)) + @test isnothing(boundary) + @test all(iszero, system.cache.surface_normal) + + system, boundary, semi, + ode = create_fluid_system(coordinates, velocity, mass, density, particle_spacing, + nothing; + smoothing_length=0.15, + surface_normal_method=ColorfieldSurfaceNormal(interface_threshold=0.0)) + @test any(!iszero, system.cache.surface_normal) + + system, boundary, semi, + ode = create_fluid_system(coordinates, velocity, mass, density, particle_spacing, + nothing; + smoothing_length=0.15, + surface_normal_method=ColorfieldSurfaceNormal(interface_threshold=0.0, + ideal_density_threshold=0.01)) + @test isnothing(boundary) + @test all(iszero, system.cache.surface_normal) + + @test_throws ArgumentError create_fluid_system(coordinates, velocity, mass, density, + particle_spacing, nothing; + smoothing_length=0.15, color_value=0, + wall=true, walldistance=particle_spacing) +end + function compute_and_test_surface_values(system, semi, ode; NDIMS=2) v0_ode, u0_ode = ode.u0.x v = TrixiParticles.wrap_v(v0_ode, system, semi) diff --git a/test/schemes/fluid/surface_tension.jl b/test/schemes/fluid/surface_tension.jl index 687de166ef..d83414112f 100644 --- a/test/schemes/fluid/surface_tension.jl +++ b/test/schemes/fluid/surface_tension.jl @@ -28,18 +28,17 @@ ideal_density_threshold=0.0f0) isa ColorfieldSurfaceNormal{Float32} - invalid_thresholds = ((NaN, 0.01, 0.0), - (-Inf, 0.01, 0.0), - (0.1, Inf, 0.0), - (0.1, 0.01, 1.0im), - ("invalid", 0.01, 0.0)) - for (boundary_threshold, interface_threshold, density_threshold) in - invalid_thresholds - - @test_throws ArgumentError ColorfieldSurfaceNormal(; - boundary_contact_threshold=boundary_threshold, - interface_threshold, - ideal_density_threshold=density_threshold) + for normal_method in + (() -> ColorfieldSurfaceNormal(boundary_contact_threshold=-0.1), + () -> ColorfieldSurfaceNormal(boundary_contact_threshold=1.1), + () -> ColorfieldSurfaceNormal(boundary_contact_threshold=NaN), + () -> ColorfieldSurfaceNormal(boundary_contact_threshold="invalid"), + () -> ColorfieldSurfaceNormal(interface_threshold=-0.1), + () -> ColorfieldSurfaceNormal(interface_threshold=Inf), + () -> ColorfieldSurfaceNormal(ideal_density_threshold=-0.1), + () -> ColorfieldSurfaceNormal(ideal_density_threshold=1.1), + () -> ColorfieldSurfaceNormal(interface_threshold=1.0im)) + @test_throws ArgumentError normal_method() end end @@ -58,7 +57,7 @@ state_equation=StateEquationCole(sound_speed=10.0, reference_density=1.0, exponent=1), - surface_tension, color_value=0) + surface_tension, color_value=1) edac = EntropicallyDampedSPHSystem(initial_condition; smoothing_kernel, smoothing_length, sound_speed=10.0, density_calculator=SummationDensity(), @@ -107,7 +106,8 @@ color_value=0) @test full_akinci.surface_normal_method isa ColorfieldSurfaceNormal @test haskey(full_akinci.cache, :surface_normal) - @test_throws ArgumentError TrixiParticles.check_system_color((full_akinci, wcsph)) + @test isnothing(TrixiParticles.check_system_color((full_akinci, wcsph))) + @test_throws ArgumentError TrixiParticles.check_system_color((full_akinci, edac)) end @testset "zero Morris coefficient does not restrict the time step" begin @@ -291,15 +291,25 @@ density_calculator=density_calc, state_equation=eq_state, surface_tension=SurfaceTensionMomentumMorris(surface_tension_coefficient=1.0), - surface_normal_method=ColorfieldSurfaceNormal(interface_threshold=0.1, - ideal_density_threshold=0.9), - reference_particle_spacing=1.0,) + surface_normal_method=ColorfieldSurfaceNormal(interface_threshold=0.0), + reference_particle_spacing=1.0) # 4. Verify Cache Contains Necessary Fields @test haskey(system.cache, :delta_s) @test haskey(system.cache, :surface_normal) @test haskey(system.cache, :stress_tensor) + # Filtering retains the raw gradient magnitude as the surface delta before + # normalizing the direction used in the stress tensor. + system.cache.surface_normal .= [3.0 0.0; + 4.0 2.0] + system.cache.neighbor_count .= 10 + TrixiParticles.remove_invalid_normals!(system, system.surface_tension, + system.surface_normal_method) + @test system.cache.surface_normal[:, 1] ≈ [0.6, 0.8] + @test system.cache.surface_normal[:, 2] ≈ [0.0, 1.0] + @test system.cache.delta_s ≈ [5.0, 2.0] + # 5. Manually Populate `delta_s` and `surface_normal` system.cache.delta_s .= [1.0, 2.0] system.cache.surface_normal .= hcat([1.0, 0.0], [1 / sqrt(2), 1 / sqrt(2)]) diff --git a/test/systems/iisph_system.jl b/test/systems/iisph_system.jl index 2bcd740552..11e026979d 100644 --- a/test/systems/iisph_system.jl +++ b/test/systems/iisph_system.jl @@ -41,6 +41,7 @@ min_iterations = min_iterations_[i] max_iterations = max_iterations_[i] time_step = time_steps_[i] + color_value = i TrixiParticles.ndims(::Val{:smoothing_kernel}) = i + 1 smoothing_kernel2 = Val(:smoothing_kernel2) # The wrong dimension. 2 -> 3, 3 -> 2. @@ -51,7 +52,7 @@ system = ImplicitIncompressibleSPHSystem(initial_condition; smoothing_kernel, smoothing_length, reference_density, omega, max_error, min_iterations, - max_iterations, time_step) + max_iterations, time_step, color_value) # Constructor copies input fields, applies defaults, and respects the requested dimensionality @test system isa ImplicitIncompressibleSPHSystem{NDIMS} @@ -67,6 +68,7 @@ @test system.min_iterations == min_iterations @test system.max_iterations == max_iterations @test system.time_step == time_step + @test system.cache.color == color_value @test length(system.density) == size(coordinates, 2) # A too-short acceleration vector triggers dimension validation diff --git a/test/systems/packing_system.jl b/test/systems/packing_system.jl index 20c915fb4c..3c165c40e0 100644 --- a/test/systems/packing_system.jl +++ b/test/systems/packing_system.jl @@ -83,6 +83,9 @@ smoothing_length_interpolation=1.0, background_pressure=1.0) + @test isnothing(TrixiParticles.surface_normal_method(system)) + @test !TrixiParticles.contributes_to_colorfield(system) + u = copy(initial_condition.coordinates) TrixiParticles.constrain_particles_onto_surface!(u, system, From 013717ca112af71a20b5d9ac5073b43d0615b860 Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Thu, 6 Aug 2026 12:28:42 +0200 Subject: [PATCH 08/13] Add dimensionally consistent 2D Akinci kernels --- NEWS.md | 7 + docs/src/systems/fluid.md | 68 +++++++-- src/schemes/fluid/surface_tension.jl | 63 ++++++--- src/schemes/structure/rigid_body/system.jl | 3 +- test/schemes/fluid/surface_tension.jl | 152 +++++++++++++++++++-- 5 files changed, 253 insertions(+), 40 deletions(-) diff --git a/NEWS.md b/NEWS.md index 74dcaf1ffd..dcea75f23a 100644 --- a/NEWS.md +++ b/NEWS.md @@ -6,6 +6,13 @@ used in the Julia ecosystem. Notable changes will be documented in this file for ## Version 0.5.4 +### API Changes + +- 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. + ### Important Bugfixes - Hardened surface tension model configuration by validating coefficients and surface-normal diff --git a/docs/src/systems/fluid.md b/docs/src/systems/fluid.md index 3f0490af0f..aa458dd8cc 100644 --- a/docs/src/systems/fluid.md +++ b/docs/src/systems/fluid.md @@ -522,6 +522,30 @@ 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. +!!! note "Akinci kernels in two dimensions" + Akinci et al. published the cohesion and adhesion kernels for three dimensions. In two + dimensions, TrixiParticles.jl uses an integral-matching extension: each radial 2D kernel + has the same full-space integral as its published 3D counterpart. This convention is not + part of the original model, but gives both kernels dimensions of ``L^{-d}`` in ``d`` + dimensions. Their products with particle mass are therefore independent of resolution at + a fixed smoothing-length-to-spacing ratio. Akinci surface tension is supported in two and + three dimensions only. Integral matching removes the resolution dependence of the pair + kernels, but does not turn their numerical coefficients into physical values in N/m. + + To preserve the pairwise cohesion and adhesion contributions from a previous 2D + configuration that used the 3D normalizations, scale the coefficients at its + compact-support radius ``h_c`` as + + ```math + \sigma_{\mathrm{new}} = \frac{627}{790h_c}\sigma_{\mathrm{old}}, \qquad + \beta_{\mathrm{new}} = \frac{42}{65h_c}\beta_{\mathrm{old}}. + ``` + + The migrated coefficients can then be held fixed when changing the resolution. Since + `SurfaceTensionAkinci` uses ``\sigma`` for both cohesion and the unchanged curvature term, + this migration also changes their relative weight; full-model configurations may require + additional calibration. + ### [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: @@ -539,18 +563,28 @@ It is defined by the distance between particles and the support radius ``h_c``, Mathematically: ```math -F_{\text{cohesion}} = -\sigma m_b C(r) \frac{r}{\Vert r \Vert}, +F_{\text{cohesion}} = -\sigma m_b C_d(r) \frac{r}{\Vert r \Vert}, ``` -where ``C(r)``, the cohesion kernel, is defined as: +where the dimension-dependent cohesion kernel is ```math -C(r)=\frac{32}{\pi h_c^9} +C_d(r)=\frac{K_d}{h_c^{d+6}} \begin{cases} (h_c-r)^3 r^3, & \text{if } 2r > h_c, \\ -2(h_c-r)^3 r^3 - \frac{h^6}{64}, & \text{if } r > 0 \text{ and } 2r \leq h_c, \\ +2(h_c-r)^3 r^3 - \frac{h_c^6}{64}, & \text{if } r > 0 \text{ and } 2r \leq h_c, \\ 0, & \text{otherwise.} \end{cases} +\qquad +K_2=\frac{25280}{627\pi}, \quad K_3=\frac{32}{\pi}. +``` + +The 3D constant is the published normalization. The 2D constant is chosen such that + +```math +\int_{\mathbb{R}^2} C_2(\Vert\bm{r}\Vert)\,\mathrm{d}A += \int_{\mathbb{R}^3} C_3(\Vert\bm{r}\Vert)\,\mathrm{d}V += \frac{79}{336}. ``` #### Surface area minimization force @@ -572,22 +606,36 @@ This force models the interaction between fluid and solid boundaries, simulating It uses a custom kernel with a peak at 0.75 times the support radius: ```math -F_{\text{adhesion}} = -\beta m_b A(r) \frac{r}{\Vert r \Vert}, +F_{\text{adhesion}} = -\beta m_b A_d(r) \frac{r}{\Vert r \Vert}, ``` -where ``A(r)`` is the adhesion kernel: +where the dimension-dependent adhesion kernel is ```math -A(r) = \frac{0.007}{h_c^{3.25}} +A_d(r) = \frac{b_d}{h_c^{d+1/4}} \begin{cases} \sqrt[4]{-\frac{4r^2}{h_c} + 6r - 2h_c}, & \text{if } 2r > h_c \text{ and } r \leq h_c, \\ 0, & \text{otherwise.} \end{cases} +\qquad +b_2=\frac{13}{1200}, \quad b_3=0.007. +``` + +Again, ``b_3`` is the published value and ``b_2`` matches the full-space integrals. In terms +of the dimensionless radial moments + +```math +J_d = \int_{1/2}^{1} q^{d-1}\left[2(1-q)(2q-1)\right]^{1/4}\,\mathrm{d}q, +``` + +the beta-function identities + +```math +J_2 = \frac{3}{8}B\!\left(\frac{5}{4},\frac{5}{4}\right), \qquad +J_3 = \frac{65}{224}B\!\left(\frac{5}{4},\frac{5}{4}\right) ``` -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. +give ``J_3/J_2=65/84`` and thus ``b_2=2b_3J_3/J_2=13/1200``. --- diff --git a/src/schemes/fluid/surface_tension.jl b/src/schemes/fluid/surface_tension.jl index 8ea5b0ae95..18cf7dcb8a 100644 --- a/src/schemes/fluid/surface_tension.jl +++ b/src/schemes/fluid/surface_tension.jl @@ -16,9 +16,9 @@ end This model only implements the cohesion force of the Akinci [Akinci2013](@cite) surface tension model. It does not require a surface-normal method. -The published Akinci cohesion kernel uses a three-dimensional normalization. In two-dimensional -simulations, `surface_tension_coefficient` is therefore an empirical numerical parameter and -may need to be adjusted when changing the resolution. +The three-dimensional cohesion kernel uses the normalization published by Akinci et al. In two +dimensions, TrixiParticles.jl uses an integral-matched extension that is independent of particle +resolution. See [`surface_tension`](@ref) for more details. @@ -44,9 +44,9 @@ principles outlined by Akinci [Akinci2013](@cite). This model is instrumental in behaviors of fluid surfaces, such as droplet formation and the dynamics of merging or separation, by utilizing intra-particle forces. -The published Akinci cohesion kernel uses a three-dimensional normalization. In two-dimensional -simulations, `surface_tension_coefficient` is therefore an empirical numerical parameter and -may need to be adjusted when changing the resolution. +The three-dimensional cohesion and adhesion kernels use the normalizations published by Akinci +et al. In two dimensions, TrixiParticles.jl uses integral-matched extensions that are independent +of particle resolution. See [`surface_tension`](@ref) for more details. @@ -91,6 +91,15 @@ function create_cache_surface_tension(surface_tension, ELTYPE, NDIMS, nparticles return (;) end +function create_cache_surface_tension(::AkinciTypeSurfaceTension, ELTYPE, NDIMS, + nparticles) + if NDIMS != 2 && NDIMS != 3 + throw(ArgumentError("Akinci surface tension is only supported in two and three dimensions")) + end + + return (;) +end + function create_cache_surface_tension(::SurfaceTensionMorris, ELTYPE, NDIMS, nparticles) curvature = Array{ELTYPE, 1}(undef, nparticles) return (; curvature) @@ -145,13 +154,34 @@ end # By using the `@fastpow` macro, we are consciously trading off some precision in the result # for enhanced computational speed. This is especially useful in scenarios where performance # is a higher priority than exact precision. +@inline function cohesion_kernel_normalization_akinci(support_radius, ::Val{2}) + normalization = oftype(support_radius, 25280 / (627 * pi)) + return (normalization / support_radius) / support_radius +end + +@inline function cohesion_kernel_normalization_akinci(support_radius, ::Val{3}) + normalization = oftype(support_radius, 32 / pi) + return ((normalization / support_radius) / support_radius) / support_radius +end + +@inline function adhesion_kernel_normalization_akinci(support_radius, ::Val{2}) + normalization = oftype(support_radius, 13 / 1200) + return (normalization / support_radius) / support_radius +end + +@inline function adhesion_kernel_normalization_akinci(support_radius, ::Val{3}) + normalization = oftype(support_radius, 0.007) + return ((normalization / support_radius) / support_radius) / support_radius +end + @fastpow @inline function cohesion_force_akinci(surface_tension, support_radius, m_b, - pos_diff, distance) + pos_diff, distance, dimensions) (; surface_tension_coefficient) = surface_tension distance >= support_radius && return zero(pos_diff) - # Eq. 2 in dimensionless form avoids scale-dependent powers up to `support_radius^9`. + # Eq. 2 in dimensionless form, with the published normalization in 3D and an + # integral-matched one in 2D. # We only reach this function when `sqrt(eps()) < distance <= support_radius` normalized_distance = distance / support_radius if normalized_distance > one(normalized_distance) / 2 @@ -163,10 +193,7 @@ end C = 2 * (1 - normalized_distance)^3 * normalized_distance^3 - one(normalized_distance) / 64 end - normalization = oftype(support_radius, 32 / pi) - normalization = ((normalization / support_radius) / support_radius) / - support_radius - C *= normalization + C *= cohesion_kernel_normalization_akinci(support_radius, dimensions) # Eq. 1 in acceleration form cohesion_force = -surface_tension_coefficient * m_b * C * pos_diff / distance @@ -175,7 +202,7 @@ end end @inline function adhesion_force_akinci(surface_tension, support_radius, m_b, pos_diff, - distance, adhesion_coefficient) + distance, adhesion_coefficient, dimensions) distance >= support_radius && return zero(pos_diff) distance <= 0.5 * support_radius && return zero(pos_diff) @@ -184,8 +211,7 @@ end normalized_distance = distance / support_radius radicand = 2 * (2 * normalized_distance - 1) * (1 - normalized_distance) fourth_root = sqrt(sqrt(max(zero(radicand), radicand))) - normalization = convert(typeof(support_radius), 0.007) - normalization = ((normalization / support_radius) / support_radius) / support_radius + normalization = adhesion_kernel_normalization_akinci(support_radius, dimensions) A = normalization * fourth_root # Eq. 6 in acceleration form with `m_b` being the boundary mass calculated as @@ -222,7 +248,7 @@ end dv_particle[] += surface_tension_correction * cohesion_force_akinci(surface_tension_a, support_radius, m_b, - pos_diff, distance) + pos_diff, distance, Val(ndims(particle_system))) return dv_particle end @@ -249,7 +275,7 @@ end dv_particle[] += surface_tension_correction * cohesion_force_akinci(surface_tension_a, support_radius, m_b, - pos_diff, distance) + pos_diff, distance, Val(ndims(particle_system))) dv_particle[] -= surface_tension_correction * surface_tension_coefficient * (n_a - n_b) * smoothing_length_ @@ -366,7 +392,8 @@ end support_radius = compact_support(particle_system.smoothing_kernel, smoothing_length(particle_system, particle)) dv_particle[] += adhesion_force_akinci(surface_tension, support_radius, m_b, pos_diff, - distance, adhesion_coefficient) + distance, adhesion_coefficient, + Val(ndims(particle_system))) return dv_particle end diff --git a/src/schemes/structure/rigid_body/system.jl b/src/schemes/structure/rigid_body/system.jl index 4f3b79090e..9b9ff16094 100644 --- a/src/schemes/structure/rigid_body/system.jl +++ b/src/schemes/structure/rigid_body/system.jl @@ -302,7 +302,8 @@ end smoothing_length(particle_system, particle)) dv_particle[] += adhesion_force_akinci(surface_tension, support_radius, m_b, - pos_diff, distance, adhesion_coefficient) + pos_diff, distance, adhesion_coefficient, + Val(ndims(particle_system))) return dv_particle end diff --git a/test/schemes/fluid/surface_tension.jl b/test/schemes/fluid/surface_tension.jl index d83414112f..4db98816e4 100644 --- a/test/schemes/fluid/surface_tension.jl +++ b/test/schemes/fluid/surface_tension.jl @@ -150,35 +150,40 @@ # Additional digits have been accepted from the actual calculation. test_distance = 0.1 val = TrixiParticles.cohesion_force_akinci(surface_tension, support_radius, m_b, - pos_diff, test_distance) * test_distance + pos_diff, test_distance, Val(3)) * + test_distance @test isapprox(val[1], 0.1443038770421044, atol=6e-15) @test isapprox(val[2], 0.1443038770421044, atol=6e-15) # Maximum repulsion force test_distance = 0.01 max = TrixiParticles.cohesion_force_akinci(surface_tension, support_radius, m_b, - pos_diff, test_distance) * test_distance + pos_diff, test_distance, Val(3)) * + test_distance @test isapprox(max[1], 0.15913517632298307, atol=6e-15) @test isapprox(max[2], 0.15913517632298307, atol=6e-15) # Near 0 test_distance = 0.2725 zero = TrixiParticles.cohesion_force_akinci(surface_tension, support_radius, m_b, - pos_diff, test_distance) * test_distance + pos_diff, test_distance, Val(3)) * + test_distance @test isapprox(zero[1], 0.0004360543645195717, atol=6e-15) @test isapprox(zero[2], 0.0004360543645195717, atol=6e-15) # Maximum attraction force test_distance = 0.5 maxa = TrixiParticles.cohesion_force_akinci(surface_tension, support_radius, m_b, - pos_diff, test_distance) * test_distance + pos_diff, test_distance, Val(3)) * + test_distance @test isapprox(maxa[1], -0.15915494309189535, atol=6e-15) @test isapprox(maxa[2], -0.15915494309189535, atol=6e-15) # Should be 0 test_distance = 1.0 zero = TrixiParticles.cohesion_force_akinci(surface_tension, support_radius, m_b, - pos_diff, test_distance) * test_distance + pos_diff, test_distance, Val(3)) * + test_distance @test isapprox(zero[1], 0.0, atol=6e-15) @test isapprox(zero[2], 0.0, atol=6e-15) @@ -207,14 +212,14 @@ # Additional digits have been accepted from the actual calculation. test_distance = 0.1 zero = TrixiParticles.adhesion_force_akinci(surface_tension, support_radius, m_b, - pos_diff, test_distance, 1.0) * + pos_diff, test_distance, 1.0, Val(3)) * test_distance @test isapprox(zero[1], 0.0, atol=6e-15) @test isapprox(zero[2], 0.0, atol=6e-15) test_distance = 0.5 zero = TrixiParticles.adhesion_force_akinci(surface_tension, support_radius, m_b, - pos_diff, test_distance, 1.0) * + pos_diff, test_distance, 1.0, Val(3)) * test_distance @test isapprox(zero[1], 0.0, atol=6e-15) @test isapprox(zero[2], 0.0, atol=6e-15) @@ -222,7 +227,7 @@ # Near 0 test_distance = 0.51 zero = TrixiParticles.adhesion_force_akinci(surface_tension, support_radius, m_b, - pos_diff, test_distance, 1.0) * + pos_diff, test_distance, 1.0, Val(3)) * test_distance @test isapprox(zero[1], -0.002619160170741761, atol=6e-15) @test isapprox(zero[2], -0.002619160170741761, atol=6e-15) @@ -230,7 +235,7 @@ # Maximum adhesion force test_distance = 0.75 max = TrixiParticles.adhesion_force_akinci(surface_tension, support_radius, m_b, - pos_diff, test_distance, 1.0) * + pos_diff, test_distance, 1.0, Val(3)) * test_distance @test isapprox(max[1], -0.004949747468305833, atol=6e-15) @test isapprox(max[2], -0.004949747468305833, atol=6e-15) @@ -238,7 +243,7 @@ # Should be 0 test_distance = 1.0 zero = TrixiParticles.adhesion_force_akinci(surface_tension, support_radius, m_b, - pos_diff, test_distance, 1.0) * + pos_diff, test_distance, 1.0, Val(3)) * test_distance @test isapprox(zero[1], 0.0, atol=6e-15) @test isapprox(zero[2], 0.0, atol=6e-15) @@ -248,7 +253,7 @@ near_support = TrixiParticles.adhesion_force_akinci(surface_tension, support_radius_f32, 1.0f0, Float32[1, 0], distance_f32, - 1.0f0) + 1.0f0, Val(3)) @test eltype(near_support) == Float32 @test all(isfinite, near_support) @test 0 < norm(near_support) < eps(Float32) @@ -266,6 +271,131 @@ end end + @testset "two-dimensional Akinci kernels" begin + surface_tension = SurfaceTensionAkinci(surface_tension_coefficient=1.0) + support_radius = 1.0 + cohesion_normalization = 25280 / (627 * pi) + + for distance in (0.25, 0.75) + pos_diff = SVector(distance, 0.0) + shape = if distance > 0.5 * support_radius + (support_radius - distance)^3 * distance^3 + else + 2 * (support_radius - distance)^3 * distance^3 - support_radius^6 / 64 + end + expected = -cohesion_normalization * shape * pos_diff / distance + force = TrixiParticles.cohesion_force_akinci(surface_tension, support_radius, + 1.0, pos_diff, distance, Val(2)) + @test isapprox(force, expected; rtol=5eps(), atol=5eps()) + end + + distance = 0.75 + pos_diff = SVector(distance, 0.0) + radicand = -4 * distance^2 / support_radius + 6 * distance - + 2 * support_radius + expected = -(13 / 1200) * radicand^(1 / 4) * pos_diff / distance + force = TrixiParticles.adhesion_force_akinci(surface_tension, support_radius, 1.0, + pos_diff, distance, 1.0, Val(2)) + @test isapprox(force, expected; rtol=5eps(), atol=5eps()) + + surface_tension_f32 = SurfaceTensionAkinci(surface_tension_coefficient=1.0f0) + distance_f32 = 0.75f0 + pos_diff_f32 = SVector(distance_f32, 0.0f0) + cohesion_f32 = TrixiParticles.cohesion_force_akinci(surface_tension_f32, 1.0f0, + 1.0f0, pos_diff_f32, + distance_f32, Val(2)) + adhesion_f32 = TrixiParticles.adhesion_force_akinci(surface_tension_f32, 1.0f0, + 1.0f0, pos_diff_f32, + distance_f32, 1.0f0, Val(2)) + @test eltype(cohesion_f32) == Float32 + @test eltype(adhesion_f32) == Float32 + @test all(isfinite, cohesion_f32) + @test all(isfinite, adhesion_f32) + + for dimensions in (1, 4) + @test_throws ArgumentError TrixiParticles.create_cache_surface_tension(surface_tension, + Float64, + dimensions, + 1) + end + end + + @testset "Akinci kernel resolution scaling" begin + surface_tension = SurfaceTensionAkinci(surface_tension_coefficient=0.8) + adhesion_coefficient = 0.6 + + function forces(scale, dimensions::Val{NDIMS}) where {NDIMS} + support_radius = scale + distance = 0.75 * support_radius + pos_diff = SVector{NDIMS}(ntuple(i -> i == 1 ? distance : zero(distance), + NDIMS)) + mass = scale^NDIMS + cohesion = TrixiParticles.cohesion_force_akinci(surface_tension, + support_radius, mass, + pos_diff, distance, dimensions) + adhesion = TrixiParticles.adhesion_force_akinci(surface_tension, + support_radius, mass, + pos_diff, distance, + adhesion_coefficient, + dimensions) + return cohesion, adhesion + end + + for dimensions in (Val(2), Val(3)) + reference_cohesion, reference_adhesion = forces(1.0, dimensions) + for scale in (0.25, 0.5, 2.0, 4.0) + cohesion, adhesion = forces(scale, dimensions) + @test isapprox(cohesion, reference_cohesion; rtol=5eps(), atol=5eps()) + @test isapprox(adhesion, reference_adhesion; rtol=5eps(), atol=5eps()) + end + end + end + + @testset "Akinci kernel integral matching" begin + surface_tension = SurfaceTensionAkinci(surface_tension_coefficient=1.0) + support_radius = 1.3 + + function pos_diff_at_radius(radius, ::Val{NDIMS}) where {NDIMS} + return SVector{NDIMS}(ntuple(i -> i == 1 ? radius : zero(radius), NDIMS)) + end + + function integrate_cohesion(dimensions::Val{NDIMS}) where {NDIMS} + radial_integral, + _ = quadgk(0.0, support_radius / 2, support_radius; + rtol=1e-13) do radius + pos_diff = pos_diff_at_radius(radius, dimensions) + force = TrixiParticles.cohesion_force_akinci(surface_tension, + support_radius, 1.0, + pos_diff, radius, dimensions) + return radius^(NDIMS - 1) * -force[1] + end + surface_measure = NDIMS == 2 ? 2pi : 4pi + return surface_measure * radial_integral + end + + function integrate_adhesion(dimensions::Val{NDIMS}) where {NDIMS} + radial_integral, + _ = quadgk(support_radius / 2, support_radius; + rtol=1e-13) do radius + pos_diff = pos_diff_at_radius(radius, dimensions) + force = TrixiParticles.adhesion_force_akinci(surface_tension, + support_radius, 1.0, + pos_diff, radius, 1.0, + dimensions) + return radius^(NDIMS - 1) * -force[1] + end + surface_measure = NDIMS == 2 ? 2pi : 4pi + return surface_measure * radial_integral + end + + cohesion_2d = integrate_cohesion(Val(2)) + cohesion_3d = integrate_cohesion(Val(3)) + @test isapprox(cohesion_2d, 79 / 336; rtol=1e-12) + @test isapprox(cohesion_3d, 79 / 336; rtol=1e-12) + @test isapprox(integrate_adhesion(Val(2)), integrate_adhesion(Val(3)); + rtol=1e-12) + end + @testset "compute_stress_tensors! (MomentumMorris)" begin # 1. Define Minimal Initial Condition with 2 Particles in 2D coords = [0.0 1.0; From 07ddec519e0b0f35897ac4a7ddd29b5b7424c631 Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Fri, 14 Aug 2026 16:49:05 +0200 Subject: [PATCH 09/13] Update Float32 tests for dimensional kernels --- test/schemes/fluid/surface_tension.jl | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/test/schemes/fluid/surface_tension.jl b/test/schemes/fluid/surface_tension.jl index 4db98816e4..6fdc03b5f5 100644 --- a/test/schemes/fluid/surface_tension.jl +++ b/test/schemes/fluid/surface_tension.jl @@ -192,7 +192,8 @@ distance = 0.75f0 * support_radius force = TrixiParticles.cohesion_force_akinci(surface_tension_f32, support_radius, 1.0f0, - Float32[distance, 0], distance) + Float32[distance, 0], distance, + Val(3)) expected = Float32(-32 / pi * (1 - 0.75)^3 * 0.75^3 / Float64(support_radius)^3) @test eltype(force) == Float32 @@ -263,7 +264,7 @@ force = TrixiParticles.adhesion_force_akinci(surface_tension, support_radius, 1.0f0, Float32[distance, 0], distance, - 1.0f0) + 1.0f0, Val(3)) expected = Float32(-0.007 / Float64(support_radius)^3 / sqrt(2)) @test all(isfinite, force) @test isapprox(force[1], expected; rtol=4eps(Float32)) From 56d4873c1ac78d8337bdadd143209853302950cb Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Thu, 6 Aug 2026 14:05:10 +0200 Subject: [PATCH 10/13] Fix Akinci free-surface correction with continuity density --- NEWS.md | 6 +- docs/src/systems/fluid.md | 29 ++++ src/general/corrections.jl | 26 +++- .../fluid/weakly_compressible_sph/rhs.jl | 20 ++- .../fluid/weakly_compressible_sph/system.jl | 16 ++ test/schemes/fluid/surface_tension.jl | 142 ++++++++++++++++++ 6 files changed, 229 insertions(+), 10 deletions(-) diff --git a/NEWS.md b/NEWS.md index dcea75f23a..9f2abd14cd 100644 --- a/NEWS.md +++ b/NEWS.md @@ -13,13 +13,17 @@ used in the Julia ecosystem. Notable changes will be documented in this file for 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. -### Important Bugfixes +### 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 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 diff --git a/docs/src/systems/fluid.md b/docs/src/systems/fluid.md index aa458dd8cc..3bb6c4ece9 100644 --- a/docs/src/systems/fluid.md +++ b/docs/src/systems/fluid.md @@ -600,6 +600,35 @@ where ``n_a`` and ``n_b`` are the unnormalized color gradients of the interactin ``h_a`` is the smoothing length of particle ``a``. The factor ``h_a`` makes the color-normal term dimensionless, consistent with the Akinci formulation. +#### Combined-force correction + +To compensate for particle-neighborhood deficiency at a free surface, the cohesion and +curvature contributions are multiplied by the symmetric factor + +```math +K_{ab} = \frac{2\rho_0}{\rho_a + \rho_b}. +``` + +[`AkinciFreeSurfaceCorrection`](@ref) implements this factor for the combined fluid-fluid +surface tension force. Section 4 of [Akinci et al. (2013)](@cite Akinci2013) also applies the +factor to viscosity for the same particle-deficiency reason. It does not modify pressure or wall +adhesion forces. + +The published correction assumes that the density estimate reflects missing neighbors. With +[`SummationDensity`](@ref), ``\rho_a`` and ``\rho_b`` in ``K_{ab}`` are the current densities. +For [`ContinuityDensity`](@ref) in a [`WeaklyCompressibleSPHSystem`](@ref), TrixiParticles.jl +reconstructs the auxiliary densities + +```math +\widetilde{\rho}_a = \sum_b m_b W_{ab} +``` + +and uses ``\widetilde{\rho}_a`` and ``\widetilde{\rho}_b`` only in ``K_{ab}``. Pressure and all +other density-dependent terms continue to use the integrated continuity density. The auxiliary +sum includes dummy boundary particles, so a wall that completes the particle neighborhood is not +misclassified as a free surface. This extension makes the correction independent of the selected +density calculator at the cost of one additional density-summation neighbor loop per update stage. + #### Wall adhesion force This force models the interaction between fluid and solid boundaries, simulating adhesion effects at walls. diff --git a/src/general/corrections.jl b/src/general/corrections.jl index d97eeaf341..2b2d72a06e 100644 --- a/src/general/corrections.jl +++ b/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,13 @@ 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), TrixiParticles.jl +reconstructs an auxiliary summation density for this correction 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 +46,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 +464,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 diff --git a/src/schemes/fluid/weakly_compressible_sph/rhs.jl b/src/schemes/fluid/weakly_compressible_sph/rhs.jl index 8ceb29d27f..93f9224843 100644 --- a/src/schemes/fluid/weakly_compressible_sph/rhs.jl +++ b/src/schemes/fluid/weakly_compressible_sph/rhs.jl @@ -71,10 +71,15 @@ function interact!(dv, v_particle_system, u_particle_system, # Determine correction factors. # This can usually be ignored, as these are all 1 when no correction is used. + 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, - rho_a, rho_b) + correction_rho_a, + correction_rho_b) # For `ContinuityDensity` without correction, this is equivalent to # dv_pressure = -m_b * (p_a + p_b) / (rho_a * rho_b) * grad_kernel @@ -126,6 +131,19 @@ function interact!(dv, v_particle_system, u_particle_system, return dv end +@inline function correction_density(::AkinciFreeSurfaceCorrection, + system::WeaklyCompressibleSPHSystem, + particle, density) + if system.density_calculator isa ContinuityDensity && + haskey(system.cache, :kernel_summation_density) + return @inbounds system.cache.kernel_summation_density[particle] + end + + return density +end + +@inline correction_density(correction, system, particle, density) = density + @propagate_inbounds function neighbor_pressure(v_neighbor_system, neighbor_system, neighbor, p_a) return current_pressure(v_neighbor_system, neighbor_system, neighbor) diff --git a/src/schemes/fluid/weakly_compressible_sph/system.jl b/src/schemes/fluid/weakly_compressible_sph/system.jl index 76cdb9893e..dafed7bdb0 100644 --- a/src/schemes/fluid/weakly_compressible_sph/system.jl +++ b/src/schemes/fluid/weakly_compressible_sph/system.jl @@ -340,6 +340,22 @@ function update_pressure!(system::WeaklyCompressibleSPHSystem, v, u, v_ode, u_od return system end +function compute_correction_values!(system::WeaklyCompressibleSPHSystem, + ::AkinciFreeSurfaceCorrection, u, + v_ode, u_ode, semi) + compute_akinci_correction_density!(system, system.density_calculator, u, u_ode, semi) + return system +end + +function compute_akinci_correction_density!(system, ::ContinuityDensity, u, u_ode, semi) + summation_density!(system, semi, u, u_ode, system.cache.kernel_summation_density) + return system +end + +function compute_akinci_correction_density!(system, ::SummationDensity, u, u_ode, semi) + return system +end + function update_final!(system::WeaklyCompressibleSPHSystem, v, u, v_ode, u_ode, semi, t; kwargs...) (; surface_tension) = system diff --git a/test/schemes/fluid/surface_tension.jl b/test/schemes/fluid/surface_tension.jl index 6fdc03b5f5..d5590d6bea 100644 --- a/test/schemes/fluid/surface_tension.jl +++ b/test/schemes/fluid/surface_tension.jl @@ -397,6 +397,148 @@ rtol=1e-12) end + @testset "Akinci free-surface correction" begin + correction = AkinciFreeSurfaceCorrection(1000.0) + @test TrixiParticles.free_surface_correction(correction, nothing, 1000.0, + 1000.0) == (1.0, 1, 1.0) + expected = 1000.0 / ((500.0 + 1000.0) / 2) + viscosity, pressure, + surface_tension = TrixiParticles.free_surface_correction(correction, nothing, + 500.0, 1000.0) + @test viscosity == expected + @test pressure == 1 + @test surface_tension == expected + @test TrixiParticles.free_surface_correction(correction, nothing, 1000.0, + 500.0) == (expected, 1, expected) + end + + @testset "Akinci ContinuityDensity reconstruction" begin + particle_spacing = 1.0 + rho0 = 1000.0 + smoothing_kernel = SchoenbergCubicSplineKernel{2}() + state_equation = StateEquationCole(sound_speed=10.0, reference_density=rho0, + exponent=1) + correction = AkinciFreeSurfaceCorrection(rho0) + fluid = RectangularShape(particle_spacing, (7, 7), (0.0, 0.0); density=rho0) + + function correction_density_values(density_calculator) + system = WeaklyCompressibleSPHSystem(fluid; smoothing_kernel, + smoothing_length=particle_spacing, + density_calculator, state_equation, + correction) + semi = Semidiscretization(system) + ode = semidiscretize(semi, (0.0, 0.01)) + v_ode, u_ode = ode.u0.x + TrixiParticles.update_systems_and_nhs(v_ode, u_ode, semi, 0.0) + density = GC.@preserve v_ode begin + v = TrixiParticles.wrap_v(v_ode, system, semi) + collect(TrixiParticles.current_density(v, system)) + end + correction_density = [TrixiParticles.correction_density(correction, system, + particle, + density[particle]) + for particle in TrixiParticles.eachparticle(system)] + return system, density, correction_density + end + + continuity_system, continuity_density, + continuity_correction_density = correction_density_values(ContinuityDensity()) + _, summation_density, + summation_correction_density = correction_density_values(SummationDensity()) + + @test all(==(rho0), continuity_density) + @test isapprox(continuity_system.cache.kernel_summation_density, + summation_density; rtol=2eps()) + @test isapprox(continuity_correction_density, + summation_correction_density; rtol=2eps()) + + coordinates = fluid.coordinates + particle_at(position) = findfirst(particle -> coordinates[:, particle] == position, + axes(coordinates, 2)) + center = particle_at([3.5, 3.5]) + face = particle_at([3.5, 0.5]) + corner = particle_at([0.5, 0.5]) + k = rho0 ./ continuity_correction_density + + @test isapprox(k[center], 1; atol=0.002) + @test k[face] > 1.15 + @test k[corner] > k[face] + + # Dummy boundary masses complete the kernel sum at a wall, so wall particles are + # not mistaken for a free surface by the reconstructed density. + tank = RectangularTank(particle_spacing, (7.0, 5.0), (7.0, 8.0), rho0; + n_layers=2, faces=(false, false, true, false)) + wall_system = WeaklyCompressibleSPHSystem(tank.fluid; smoothing_kernel, + smoothing_length=particle_spacing, + density_calculator=ContinuityDensity(), + state_equation, correction) + boundary_model = BoundaryModelDummyParticles(tank.boundary; + fluid_system=wall_system, + boundary_density_calculator=AdamiPressureExtrapolation()) + boundary_system = WallBoundarySystem(tank.boundary, boundary_model) + wall_semi = Semidiscretization(wall_system, boundary_system) + wall_ode = semidiscretize(wall_semi, (0.0, 0.01)) + TrixiParticles.update_systems_and_nhs(wall_ode.u0.x..., wall_semi, 0.0) + + wall_coordinates = tank.fluid.coordinates + wall_particle_at(position) = findfirst(particle -> wall_coordinates[:, particle] == + position, + axes(wall_coordinates, 2)) + bottom = wall_particle_at([3.5, 0.5]) + interior = wall_particle_at([3.5, 2.5]) + top = wall_particle_at([3.5, 4.5]) + reconstructed_density = wall_system.cache.kernel_summation_density + wall_k = rho0 ./ reconstructed_density + + @test isapprox(wall_k[bottom], wall_k[interior]; rtol=2eps()) + @test isapprox(wall_k[interior], 1; atol=0.002) + @test wall_k[top] > 1.15 + end + + @testset "Akinci correction force assembly" begin + rho0 = 1000.0 + particle_spacing = 0.5 + coordinates = [0.0 0.75; 0.0 0.0] + initial_condition = InitialCondition(; coordinates, velocity=zeros(2, 2), + mass=fill(rho0 * particle_spacing^2, 2), + density=fill(rho0, 2), particle_spacing) + smoothing_kernel = WendlandC2Kernel{2}() + state_equation = StateEquationCole(sound_speed=10.0, reference_density=rho0, + exponent=1) + surface_tension = CohesionForceAkinci(surface_tension_coefficient=0.2) + + function initial_acceleration(correction) + system = WeaklyCompressibleSPHSystem(initial_condition; smoothing_kernel, + smoothing_length=0.5, + density_calculator=ContinuityDensity(), + state_equation, surface_tension, + correction, + reference_particle_spacing=particle_spacing) + semi = Semidiscretization(system) + ode = semidiscretize(semi, (0.0, 0.01)) + v_ode, u_ode = ode.u0.x + TrixiParticles.update_systems_and_nhs(v_ode, u_ode, semi, 0.0) + dv = GC.@preserve v_ode u_ode begin + v = TrixiParticles.wrap_v(v_ode, system, semi) + u = TrixiParticles.wrap_u(u_ode, system, semi) + dv_inner = zeros(eltype(v), size(v)) + TrixiParticles.interact!(dv_inner, v, u, v, u, system, system, semi) + dv_inner + end + return system, dv[1:2, :] + end + + corrected_system, + corrected_acceleration = initial_acceleration(AkinciFreeSurfaceCorrection(rho0)) + _, uncorrected_acceleration = initial_acceleration(nothing) + correction_factor = rho0 / corrected_system.cache.kernel_summation_density[1] + + @test correction_factor > 1 + @test maximum(abs, uncorrected_acceleration) > 0 + @test isapprox(corrected_acceleration, + correction_factor * uncorrected_acceleration; rtol=2eps()) + end + @testset "compute_stress_tensors! (MomentumMorris)" begin # 1. Define Minimal Initial Condition with 2 Particles in 2D coords = [0.0 1.0; From 67980d207eeff8b91ecbbf9b986933727922cce1 Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Sat, 8 Aug 2026 16:38:55 +0200 Subject: [PATCH 11/13] Complete Akinci correction support for EDAC --- NEWS.md | 11 ++- docs/src/systems/entropically_damped_sph.md | 5 + docs/src/systems/fluid.md | 16 +-- src/general/corrections.jl | 10 +- .../fluid/entropically_damped_sph/rhs.jl | 15 ++- .../fluid/entropically_damped_sph/system.jl | 25 ++++- src/schemes/fluid/surface_normal_sph.jl | 4 + .../fluid/weakly_compressible_sph/rhs.jl | 3 +- .../fluid/weakly_compressible_sph/system.jl | 17 ++++ test/schemes/fluid/surface_tension.jl | 97 ++++++++++++++++++- 10 files changed, 180 insertions(+), 23 deletions(-) diff --git a/NEWS.md b/NEWS.md index 9f2abd14cd..33d90c8902 100644 --- a/NEWS.md +++ b/NEWS.md @@ -13,6 +13,11 @@ used in the Julia ecosystem. Notable changes will be documented in this file for 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. + ### Important Bugfixes - Hardened surface tension model configuration by validating coefficients and surface-normal @@ -21,9 +26,9 @@ used in the Julia ecosystem. Notable changes will be documented in this file for - 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 without changing the integrated - pressure density. Dummy boundary masses are included so wall particles are not misclassified - as a free surface. + 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 diff --git a/docs/src/systems/entropically_damped_sph.md b/docs/src/systems/entropically_damped_sph.md index 96acbad352..0df6b23f8e 100644 --- a/docs/src/systems/entropically_damped_sph.md +++ b/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`` diff --git a/docs/src/systems/fluid.md b/docs/src/systems/fluid.md index 3bb6c4ece9..ae65c75484 100644 --- a/docs/src/systems/fluid.md +++ b/docs/src/systems/fluid.md @@ -616,18 +616,22 @@ adhesion forces. The published correction assumes that the density estimate reflects missing neighbors. With [`SummationDensity`](@ref), ``\rho_a`` and ``\rho_b`` in ``K_{ab}`` are the current densities. -For [`ContinuityDensity`](@ref) in a [`WeaklyCompressibleSPHSystem`](@ref), TrixiParticles.jl +For [`ContinuityDensity`](@ref) in a [`WeaklyCompressibleSPHSystem`](@ref) or +[`EntropicallyDampedSPHSystem`](@ref), TrixiParticles.jl reconstructs the auxiliary densities ```math \widetilde{\rho}_a = \sum_b m_b W_{ab} ``` -and uses ``\widetilde{\rho}_a`` and ``\widetilde{\rho}_b`` only in ``K_{ab}``. Pressure and all -other density-dependent terms continue to use the integrated continuity density. The auxiliary -sum includes dummy boundary particles, so a wall that completes the particle neighborhood is not -misclassified as a free surface. This extension makes the correction independent of the selected -density calculator at the cost of one additional density-summation neighbor loop per update stage. +and uses ``\widetilde{\rho}_a`` and ``\widetilde{\rho}_b`` in ``K_{ab}``. For +[`SurfaceTensionAkinci`](@ref), the same reconstructed densities are used in the particle-volume +weights ``m_b / \rho_b`` of the color-field normals, matching the density-summation formulation +of the published model. Pressure and all other density-dependent terms continue to use the +integrated continuity density. The auxiliary sum includes dummy boundary particles, so a wall +that completes the particle neighborhood is not misclassified as a free surface. This extension +makes the correction independent of the selected density calculator at the cost of one additional +density-summation neighbor loop per update stage. #### Wall adhesion force diff --git a/src/general/corrections.jl b/src/general/corrections.jl index 2b2d72a06e..dcbd293026 100644 --- a/src/general/corrections.jl +++ b/src/general/corrections.jl @@ -21,10 +21,12 @@ this value will be about ~1.5 for particles at the free surface and can then be 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), TrixiParticles.jl -reconstructs an auxiliary summation density for this correction 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. +[`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. diff --git a/src/schemes/fluid/entropically_damped_sph/rhs.jl b/src/schemes/fluid/entropically_damped_sph/rhs.jl index 2ea1e4dc4d..99adb7352a 100644 --- a/src/schemes/fluid/entropically_damped_sph/rhs.jl +++ b/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, diff --git a/src/schemes/fluid/entropically_damped_sph/system.jl b/src/schemes/fluid/entropically_damped_sph/system.jl index 028b08090f..dc5eefe514 100644 --- a/src/schemes/fluid/entropically_damped_sph/system.jl +++ b/src/schemes/fluid/entropically_damped_sph/system.jl @@ -38,7 +38,10 @@ See [Entropically Damped Artificial Compressibility for SPH](@ref edac) for more from the local pressure (default: `true` when using shifting, `false` otherwise). - `buffer_size`: Number of buffer particles. This is needed when simulating with [`OpenBoundarySystem`](@ref). -- `correction`: Correction method used for this system. (default: no correction, see [Corrections](@ref corrections)) +- `correction`: Correction method used for this system. (default: no correction, + see [Corrections](@ref corrections)). + [`AkinciFreeSurfaceCorrection`](@ref) scales viscosity and + Akinci surface-tension forces but not pressure. - `source_terms`: Additional source terms for this system. Has to be either `nothing` (by default), or a function of `(coords, velocity, density, pressure, t)` (which are the quantities of a single particle), returning a `Tuple` @@ -254,6 +257,12 @@ end system_correction(system::EntropicallyDampedSPHSystem) = system.correction +@inline function surface_normal_density(system::EntropicallyDampedSPHSystem, particle, + density) + return surface_normal_density(system, system.surface_tension, system.correction, + system.density_calculator, particle, density) +end + @inline function current_velocity(v, system::EntropicallyDampedSPHSystem) return view(v, 1:ndims(system), :) end @@ -300,11 +309,25 @@ function update_quantities!(system::EntropicallyDampedSPHSystem, v, u, end function update_pressure!(system::EntropicallyDampedSPHSystem, v, u, v_ode, u_ode, semi, t) + compute_akinci_correction_density!(system, system.correction, + system.density_calculator, u, u_ode, semi) + compute_surface_normal!(system, system.surface_normal_method, v, u, v_ode, u_ode, semi, t) compute_surface_delta_function!(system, system.surface_tension, semi) end +function compute_akinci_correction_density!(system, correction, density_calculator, + u, u_ode, semi) + return system +end + +function compute_akinci_correction_density!(system, ::AkinciFreeSurfaceCorrection, + density_calculator, u, u_ode, semi) + compute_akinci_correction_density!(system, density_calculator, u, u_ode, semi) + return system +end + function update_final!(system::EntropicallyDampedSPHSystem, v, u, v_ode, u_ode, semi, t; kwargs...) (; surface_tension) = system diff --git a/src/schemes/fluid/surface_normal_sph.jl b/src/schemes/fluid/surface_normal_sph.jl index ceda09aa12..764c07c8fa 100644 --- a/src/schemes/fluid/surface_normal_sph.jl +++ b/src/schemes/fluid/surface_normal_sph.jl @@ -78,6 +78,8 @@ end return extract_svector(cache.surface_normal, particle_system, particle) end +@inline surface_normal_density(system, particle, density) = density + function calc_normal!(system, neighbor_system, u_system, v, v_neighbor_system, u_neighbor_system, semi, surface_normal_method, neighbor_surface_normal_method) @@ -105,6 +107,8 @@ function calc_normal!(system::AbstractFluidSystem, neighbor_system::AbstractFlui m_b = hydrodynamic_mass(neighbor_system, neighbor) density_neighbor = current_density(v_neighbor_system, neighbor_system, neighbor) + density_neighbor = surface_normal_density(neighbor_system, neighbor, + density_neighbor) grad_kernel = smoothing_kernel_grad(system, pos_diff, distance, particle) for i in 1:ndims(system) cache.surface_normal[i, diff --git a/src/schemes/fluid/weakly_compressible_sph/rhs.jl b/src/schemes/fluid/weakly_compressible_sph/rhs.jl index 93f9224843..5cbc31fd60 100644 --- a/src/schemes/fluid/weakly_compressible_sph/rhs.jl +++ b/src/schemes/fluid/weakly_compressible_sph/rhs.jl @@ -132,7 +132,8 @@ function interact!(dv, v_particle_system, u_particle_system, end @inline function correction_density(::AkinciFreeSurfaceCorrection, - system::WeaklyCompressibleSPHSystem, + system::Union{WeaklyCompressibleSPHSystem, + EntropicallyDampedSPHSystem}, particle, density) if system.density_calculator isa ContinuityDensity && haskey(system.cache, :kernel_summation_density) diff --git a/src/schemes/fluid/weakly_compressible_sph/system.jl b/src/schemes/fluid/weakly_compressible_sph/system.jl index dafed7bdb0..fd484bdfc5 100644 --- a/src/schemes/fluid/weakly_compressible_sph/system.jl +++ b/src/schemes/fluid/weakly_compressible_sph/system.jl @@ -246,6 +246,23 @@ end system_correction(system::WeaklyCompressibleSPHSystem) = system.correction +@inline function surface_normal_density(system::WeaklyCompressibleSPHSystem, particle, + density) + return surface_normal_density(system, system.surface_tension, system.correction, + system.density_calculator, particle, density) +end + +@inline function surface_normal_density(system, surface_tension, correction, + density_calculator, particle, density) + return density +end + +@inline function surface_normal_density(system, ::SurfaceTensionAkinci, + ::AkinciFreeSurfaceCorrection, + ::ContinuityDensity, particle, density) + return @inbounds system.cache.kernel_summation_density[particle] +end + @propagate_inbounds function current_velocity(v, system::WeaklyCompressibleSPHSystem) return current_velocity(v, system.density_calculator, system) end diff --git a/test/schemes/fluid/surface_tension.jl b/test/schemes/fluid/surface_tension.jl index d5590d6bea..108aefe98c 100644 --- a/test/schemes/fluid/surface_tension.jl +++ b/test/schemes/fluid/surface_tension.jl @@ -419,13 +419,15 @@ state_equation = StateEquationCole(sound_speed=10.0, reference_density=rho0, exponent=1) correction = AkinciFreeSurfaceCorrection(rho0) + surface_tension = SurfaceTensionAkinci(surface_tension_coefficient=0.2) fluid = RectangularShape(particle_spacing, (7, 7), (0.0, 0.0); density=rho0) function correction_density_values(density_calculator) system = WeaklyCompressibleSPHSystem(fluid; smoothing_kernel, smoothing_length=particle_spacing, density_calculator, state_equation, - correction) + correction, surface_tension, + reference_particle_spacing=particle_spacing) semi = Semidiscretization(system) ode = semidiscretize(semi, (0.0, 0.01)) v_ode, u_ode = ode.u0.x @@ -443,7 +445,7 @@ continuity_system, continuity_density, continuity_correction_density = correction_density_values(ContinuityDensity()) - _, summation_density, + summation_system, summation_density, summation_correction_density = correction_density_values(SummationDensity()) @test all(==(rho0), continuity_density) @@ -451,6 +453,42 @@ summation_density; rtol=2eps()) @test isapprox(continuity_correction_density, summation_correction_density; rtol=2eps()) + @test isapprox(continuity_system.cache.surface_normal, + summation_system.cache.surface_normal; rtol=2eps()) + + function edac_correction_density_values(density_calculator) + system = EntropicallyDampedSPHSystem(fluid; smoothing_kernel, + smoothing_length=particle_spacing, + sound_speed=10.0, density_calculator, + correction, surface_tension, + reference_particle_spacing=particle_spacing) + semi = Semidiscretization(system) + ode = semidiscretize(semi, (0.0, 0.01)) + v_ode, u_ode = ode.u0.x + TrixiParticles.update_systems_and_nhs(v_ode, u_ode, semi, 0.0) + density = GC.@preserve v_ode begin + v = TrixiParticles.wrap_v(v_ode, system, semi) + collect(TrixiParticles.current_density(v, system)) + end + correction_density = [TrixiParticles.correction_density(correction, system, + particle, + density[particle]) + for particle in TrixiParticles.eachparticle(system)] + return system, density, correction_density + end + + edac_continuity_system, edac_continuity_density, + edac_continuity_correction_density = edac_correction_density_values(ContinuityDensity()) + edac_summation_system, edac_summation_density, + edac_summation_correction_density = edac_correction_density_values(SummationDensity()) + + @test all(==(rho0), edac_continuity_density) + @test isapprox(edac_continuity_system.cache.kernel_summation_density, + edac_summation_density; rtol=2eps()) + @test isapprox(edac_continuity_correction_density, + edac_summation_correction_density; rtol=2eps()) + @test isapprox(edac_continuity_system.cache.surface_normal, + edac_summation_system.cache.surface_normal; rtol=2eps()) coordinates = fluid.coordinates particle_at(position) = findfirst(particle -> coordinates[:, particle] == position, @@ -472,9 +510,11 @@ smoothing_length=particle_spacing, density_calculator=ContinuityDensity(), state_equation, correction) - boundary_model = BoundaryModelDummyParticles(tank.boundary; - fluid_system=wall_system, - boundary_density_calculator=AdamiPressureExtrapolation()) + boundary_model = BoundaryModelDummyParticles(tank.boundary.density, + tank.boundary.mass, + AdamiPressureExtrapolation(), + smoothing_kernel, particle_spacing; + state_equation, correction) boundary_system = WallBoundarySystem(tank.boundary, boundary_model) wall_semi = Semidiscretization(wall_system, boundary_system) wall_ode = semidiscretize(wall_semi, (0.0, 0.01)) @@ -539,6 +579,53 @@ correction_factor * uncorrected_acceleration; rtol=2eps()) end + @testset "EDAC Akinci correction force assembly" begin + rho0 = 1000.0 + particle_spacing = 0.5 + coordinates = [0.0 0.75; 0.0 0.0] + initial_condition = InitialCondition(; coordinates, velocity=zeros(2, 2), + mass=fill(rho0 * particle_spacing^2, 2), + density=fill(rho0, 2), particle_spacing) + smoothing_kernel = WendlandC2Kernel{2}() + surface_tension = CohesionForceAkinci(surface_tension_coefficient=0.2) + + function initial_acceleration(density_calculator, correction) + system = EntropicallyDampedSPHSystem(initial_condition; smoothing_kernel, + smoothing_length=particle_spacing, + sound_speed=10.0, density_calculator, + correction, surface_tension) + semi = Semidiscretization(system) + ode = semidiscretize(semi, (0.0, 0.01)) + v_ode, u_ode = ode.u0.x + TrixiParticles.update_systems_and_nhs(v_ode, u_ode, semi, 0.0) + + acceleration, + correction_density = GC.@preserve v_ode u_ode begin + v = TrixiParticles.wrap_v(v_ode, system, semi) + u = TrixiParticles.wrap_u(u_ode, system, semi) + dv = zeros(eltype(v), size(v)) + TrixiParticles.interact!(dv, v, u, v, u, system, system, semi) + density = TrixiParticles.current_density(v, system, 1) + dv[1:2, :], + TrixiParticles.correction_density(correction, system, 1, density) + end + return acceleration, correction_density + end + + for density_calculator in (ContinuityDensity(), SummationDensity()) + corrected_acceleration, + correction_density = initial_acceleration(density_calculator, + AkinciFreeSurfaceCorrection(rho0)) + uncorrected_acceleration, _ = initial_acceleration(density_calculator, nothing) + correction_factor = rho0 / correction_density + + @test correction_factor > 1 + @test maximum(abs, uncorrected_acceleration) > 0 + @test isapprox(corrected_acceleration, + correction_factor * uncorrected_acceleration; rtol=2eps()) + end + end + @testset "compute_stress_tensors! (MomentumMorris)" begin # 1. Define Minimal Initial Condition with 2 Particles in 2D coords = [0.0 1.0; From 17b313f4f5f8ec3c675b292762ffa332853677f9 Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Fri, 14 Aug 2026 17:20:35 +0200 Subject: [PATCH 12/13] Require free-surface correction for Akinci forces --- NEWS.md | 2 ++ docs/src/systems/fluid.md | 2 ++ examples/fluid/falling_water_spheres_2d.jl | 3 +++ .../fluid/entropically_damped_sph/system.jl | 2 ++ src/schemes/fluid/surface_tension.jl | 9 +++++++ .../fluid/weakly_compressible_sph/system.jl | 2 ++ test/schemes/fluid/surface_normal_sph.jl | 4 ++++ test/schemes/fluid/surface_tension.jl | 24 +++++++++++++++++-- test/systems/rigid_system.jl | 1 + 9 files changed, 47 insertions(+), 2 deletions(-) diff --git a/NEWS.md b/NEWS.md index 33d90c8902..204c2f73f6 100644 --- a/NEWS.md +++ b/NEWS.md @@ -8,6 +8,8 @@ used in the Julia ecosystem. Notable changes will be documented in this file for ### 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 diff --git a/docs/src/systems/fluid.md b/docs/src/systems/fluid.md index ae65c75484..4fea024520 100644 --- a/docs/src/systems/fluid.md +++ b/docs/src/systems/fluid.md @@ -521,6 +521,8 @@ the fluid-fluid surface force. Wall adhesion is controlled independently by the 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. +The full Akinci model also requires `AkinciFreeSurfaceCorrection`, which supplies the +particle-neighborhood correction in Equations 4 and 5 of Akinci et al. !!! note "Akinci kernels in two dimensions" Akinci et al. published the cohesion and adhesion kernels for three dimensions. In two diff --git a/examples/fluid/falling_water_spheres_2d.jl b/examples/fluid/falling_water_spheres_2d.jl index 1f0016ac6d..208c1b71bc 100644 --- a/examples/fluid/falling_water_spheres_2d.jl +++ b/examples/fluid/falling_water_spheres_2d.jl @@ -61,6 +61,8 @@ 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, @@ -68,6 +70,7 @@ sphere_surface_tension = EntropicallyDampedSPHSystem(sphere1; 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, diff --git a/src/schemes/fluid/entropically_damped_sph/system.jl b/src/schemes/fluid/entropically_damped_sph/system.jl index dc5eefe514..eec3526f46 100644 --- a/src/schemes/fluid/entropically_damped_sph/system.jl +++ b/src/schemes/fluid/entropically_damped_sph/system.jl @@ -124,6 +124,8 @@ function EntropicallyDampedSPHSystem(initial_condition; smoothing_kernel, smooth throw(ArgumentError("`acceleration` must be of length $NDIMS for a $(NDIMS)D problem")) end + check_akinci_correction(surface_tension, correction) + surface_normal_method = default_surface_normal_method(surface_tension, surface_normal_method) diff --git a/src/schemes/fluid/surface_tension.jl b/src/schemes/fluid/surface_tension.jl index 18cf7dcb8a..c894f9c25f 100644 --- a/src/schemes/fluid/surface_tension.jl +++ b/src/schemes/fluid/surface_tension.jl @@ -135,6 +135,15 @@ end @inline requires_surface_normal(::CohesionForceAkinci) = false @inline requires_surface_normal(::Any) = true +function check_akinci_correction(surface_tension, correction) + if surface_tension isa SurfaceTensionAkinci && + !(correction isa AkinciFreeSurfaceCorrection) + throw(ArgumentError("`SurfaceTensionAkinci` requires `AkinciFreeSurfaceCorrection`")) + end + + return surface_tension +end + function create_cache_surface_tension(::SurfaceTensionMomentumMorris, ELTYPE, NDIMS, nparticles) delta_s = Array{ELTYPE, 1}(undef, nparticles) diff --git a/src/schemes/fluid/weakly_compressible_sph/system.jl b/src/schemes/fluid/weakly_compressible_sph/system.jl index fd484bdfc5..b44f7cda59 100644 --- a/src/schemes/fluid/weakly_compressible_sph/system.jl +++ b/src/schemes/fluid/weakly_compressible_sph/system.jl @@ -132,6 +132,8 @@ function WeaklyCompressibleSPHSystem(initial_condition; smoothing_kernel, throw(ArgumentError("`ShepardKernelCorrection` cannot be used with `ContinuityDensity`")) end + check_akinci_correction(surface_tension, correction) + surface_normal_method = default_surface_normal_method(surface_tension, surface_normal_method) diff --git a/test/schemes/fluid/surface_normal_sph.jl b/test/schemes/fluid/surface_normal_sph.jl index 099387eacb..c2da842acc 100644 --- a/test/schemes/fluid/surface_normal_sph.jl +++ b/test/schemes/fluid/surface_normal_sph.jl @@ -77,6 +77,10 @@ function create_fluid_system(coordinates, velocity, mass, density, particle_spac system = WeaklyCompressibleSPHSystem(fluid; smoothing_kernel, smoothing_length, density_calculator=SummationDensity(), state_equation, + correction=surface_tension isa + SurfaceTensionAkinci ? + AkinciFreeSurfaceCorrection(first(density)) : + nothing, surface_normal_method, reference_particle_spacing=particle_spacing, surface_tension, color_value) diff --git a/test/schemes/fluid/surface_tension.jl b/test/schemes/fluid/surface_tension.jl index 108aefe98c..e4c5adb0a5 100644 --- a/test/schemes/fluid/surface_tension.jl +++ b/test/schemes/fluid/surface_tension.jl @@ -87,13 +87,32 @@ state_equation=StateEquationCole(sound_speed=10.0, reference_density=1.0, exponent=1), - surface_tension=SurfaceTensionAkinci()) + surface_tension=SurfaceTensionAkinci(), + correction=AkinciFreeSurfaceCorrection(1.0)) @test_throws ArgumentError EntropicallyDampedSPHSystem(initial_condition; smoothing_kernel, smoothing_length, sound_speed=10.0, density_calculator=SummationDensity(), - surface_tension=SurfaceTensionAkinci()) + surface_tension=SurfaceTensionAkinci(), + correction=AkinciFreeSurfaceCorrection(1.0)) + + @test_throws ArgumentError WeaklyCompressibleSPHSystem(initial_condition; + smoothing_kernel, + smoothing_length, + density_calculator=SummationDensity(), + state_equation=StateEquationCole(sound_speed=10.0, + reference_density=1.0, + exponent=1), + surface_tension=SurfaceTensionAkinci(), + reference_particle_spacing=1.0) + @test_throws ArgumentError EntropicallyDampedSPHSystem(initial_condition; + smoothing_kernel, + smoothing_length, + sound_speed=10.0, + density_calculator=SummationDensity(), + surface_tension=SurfaceTensionAkinci(), + reference_particle_spacing=1.0) full_akinci = WeaklyCompressibleSPHSystem(initial_condition; smoothing_kernel, smoothing_length, @@ -102,6 +121,7 @@ reference_density=1.0, exponent=1), surface_tension=SurfaceTensionAkinci(), + correction=AkinciFreeSurfaceCorrection(1.0), reference_particle_spacing=1.0, color_value=0) @test full_akinci.surface_normal_method isa ColorfieldSurfaceNormal diff --git a/test/systems/rigid_system.jl b/test/systems/rigid_system.jl index d549363d9d..ea0e32ca7b 100644 --- a/test/systems/rigid_system.jl +++ b/test/systems/rigid_system.jl @@ -500,6 +500,7 @@ density_calculator=SummationDensity(), state_equation, surface_tension=SurfaceTensionAkinci(surface_tension_coefficient=0.05), + correction=AkinciFreeSurfaceCorrection(fluid_density), reference_particle_spacing=particle_spacing) boundary_coordinates = reshape([1.5, 0.0], 2, 1) From 8b326b4fa89cc9e0d33283e05f7cb2db8f66809a Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Fri, 14 Aug 2026 18:44:59 +0200 Subject: [PATCH 13/13] Test Akinci forces with continuity density --- test/schemes/fluid/surface_tension.jl | 46 +++++++++++++++++++++++++++ 1 file changed, 46 insertions(+) diff --git a/test/schemes/fluid/surface_tension.jl b/test/schemes/fluid/surface_tension.jl index e4c5adb0a5..a3cbd23096 100644 --- a/test/schemes/fluid/surface_tension.jl +++ b/test/schemes/fluid/surface_tension.jl @@ -515,13 +515,59 @@ axes(coordinates, 2)) center = particle_at([3.5, 3.5]) face = particle_at([3.5, 0.5]) + face_inner = particle_at([3.5, 1.5]) corner = particle_at([0.5, 0.5]) + corner_neighbor = particle_at([1.5, 0.5]) k = rho0 ./ continuity_correction_density @test isapprox(k[center], 1; atol=0.002) @test k[face] > 1.15 @test k[corner] > k[face] + function pair_surface_acceleration(system, density, correction_density, + particle, neighbor) + pos_diff = SVector{2}(coordinates[:, particle] - coordinates[:, neighbor]) + distance = norm(pos_diff) + surface_tension_correction = TrixiParticles.free_surface_correction(correction, + system, + correction_density[particle], + correction_density[neighbor])[3] + acceleration = Ref(zero(pos_diff)) + TrixiParticles.surface_tension_force!(acceleration, surface_tension, + surface_tension, system, system, + particle, neighbor, pos_diff, distance, + density[particle], density[neighbor], + zero(pos_diff), + surface_tension_correction) + return acceleration[] + end + + # The reconstructed density must make the complete cohesion-plus-normal force agree + # with SummationDensity, not just its individual correction and normal inputs. + for (particle, neighbor) in ((face, face_inner), (corner, corner_neighbor)) + continuity_acceleration = pair_surface_acceleration(continuity_system, + continuity_density, + continuity_correction_density, + particle, neighbor) + summation_acceleration = pair_surface_acceleration(summation_system, + summation_density, + summation_correction_density, + particle, neighbor) + edac_continuity_acceleration = pair_surface_acceleration(edac_continuity_system, + edac_continuity_density, + edac_continuity_correction_density, + particle, neighbor) + edac_summation_acceleration = pair_surface_acceleration(edac_summation_system, + edac_summation_density, + edac_summation_correction_density, + particle, neighbor) + + @test !iszero(continuity_acceleration) + @test isapprox(continuity_acceleration, summation_acceleration; rtol=2eps()) + @test isapprox(edac_continuity_acceleration, edac_summation_acceleration; + rtol=2eps()) + end + # Dummy boundary masses complete the kernel sum at a wall, so wall particles are # not mistaken for a free surface by the reconstructed density. tank = RectangularTank(particle_spacing, (7.0, 5.0), (7.0, 8.0), rho0;