diff --git a/NEWS.md b/NEWS.md index ad3e44c715..80fc3dc823 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 and surface-normal + thresholds, avoiding unnecessary normal allocation for `CohesionForceAkinci`, and stabilizing + Akinci cohesion and adhesion kernels across floating-point scales. + ## Version 0.5.3 ### Features diff --git a/docs/src/systems/fluid.md b/docs/src/systems/fluid.md index c9b9860ab9..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 are essential for modeling surface tension as they provide the directionality -of forces acting at the fluid interface. They are calculated based on the particle properties and -their spatial distribution. +Surface normals characterize the local orientation of an interface. In SPH, this geometric +information can be used for interface detection and reconstruction, curvature estimation, +interfacial boundary conditions, and interfacial force models. The computed normal field is also +available for analysis and VTK output. #### Color field and gradient-based surface normals -The surface normal at a particle is derived from the color field, a scalar field assigned to particles -to distinguish between different fluid phases or between fluid and air. The color field gradients point -towards the interface, and the normalized gradient defines the surface normal direction. +The surface normal at a particle can be derived from a color field, a scalar marker used to +distinguish phases or materials. Its gradient is perpendicular to the color-field level sets and +therefore provides an interface-normal estimate; its orientation depends on the chosen color +convention. For a free surface whose exterior phase is not represented by particles, truncation of +the kernel support creates the corresponding discrete color-field gradient. -The simplest SPH formulation for a surface normal, ``n_a`` is given as +The simplest SPH approximation of an unnormalized color-field normal, ``n_a``, is ```math n_a = \sum_b m_b \frac{c_b}{\rho_b} \nabla_a W_{ab}, @@ -238,23 +241,117 @@ where: - ``\rho_b`` is the density of particle ``b``, - ``\nabla_a W_{ab}`` is the gradient of the smoothing kernel ``W_{ab}`` with respect to particle ``a``. +```@eval +using CairoMakie + +let + coordinate = range(-2.0, 2.0, length=401) + interface_width = 0.3 + colorfield = @. 0.5 * (1.0 - tanh(coordinate / interface_width)) + colorfield_gradient = @. -0.5 / interface_width / + cosh(coordinate / interface_width)^2 + + fig = Figure(size=(1000, 430), fontsize=18) + color_axis = Axis(fig[1, 1], + xlabel="signed distance s/h", ylabel="color field c", + title="Diffuse color-field transition") + gradient_axis = Axis(fig[1, 2], + xlabel="signed distance s/h", ylabel="dc/d(s/h)", + title="Color-field gradient") + + lines!(color_axis, coordinate, colorfield, color=:steelblue, linewidth=3) + lines!(gradient_axis, coordinate, colorfield_gradient, color=:darkorange, + linewidth=3) + vlines!(color_axis, [0.0], color=:black, linestyle=:dash, linewidth=2) + vlines!(gradient_axis, [0.0], color=:black, linestyle=:dash, linewidth=2) + hlines!(gradient_axis, [0.0], color=(:black, 0.35), linewidth=1) + xlims!(color_axis, extrema(coordinate)) + xlims!(gradient_axis, extrema(coordinate)) + + CairoMakie.save("colorfield_profile.png", fig) +end +``` + +![A diffuse color field and its gradient across an interface](colorfield_profile.png) + +The color field is approximately constant within either phase. Its gradient is localized in the +transition region and vanishes away from the interface. The sign of the gradient determines the +normal orientation; exchanging the two color values reverses that orientation. + #### Normalization of surface normals -The calculated normals are normalized to unit vectors: +The color-field gradient ``n_a`` is generally not a unit vector. Formulations that require only +the interface orientation use the unit normal ```math \hat{n}_a = \frac{n_a}{\Vert n_a \Vert}. ``` -Normalization ensures that the magnitude of the normals does not bias the curvature calculations or the resulting surface tension forces. +Normalization separates the interface orientation from the magnitude of the discrete color-field +gradient. In TrixiParticles, the Morris formulations use unit normals for curvature or surface-stress +calculations, while the Akinci surface-area force uses the unnormalized gradient. + +```@eval +using CairoMakie + +let + coordinate = range(-1.35, 1.35, length=241) + radius = 0.85 + interface_width = 0.1 + colorfield = [0.5 * (1.0 - tanh((hypot(x, y) - radius) / interface_width)) + for x in coordinate, y in coordinate] + + particle_spacing = 0.17 + particle_coordinates = [(x, y) for x in (-radius):particle_spacing:radius + for y in (-radius):particle_spacing:radius + if hypot(x, y) <= radius] + particle_x = first.(particle_coordinates) + particle_y = last.(particle_coordinates) + + angles = range(0.0, 2pi, length=13)[1:(end - 1)] + normal_x = -cos.(angles) + normal_y = -sin.(angles) + interface_x = radius .* cos.(angles) + interface_y = radius .* sin.(angles) + + fig = Figure(size=(760, 650), fontsize=18) + axis = Axis(fig[1, 1], aspect=DataAspect(), + xlabel="x/h", ylabel="y/h", + title="Interface orientation from the color-field gradient") + heatmap = heatmap!(axis, coordinate, coordinate, colorfield, + colormap=:viridis, colorrange=(0.0, 1.0)) + contour!(axis, coordinate, coordinate, colorfield, + levels=[0.1, 0.9], color=(:white, 0.8), linewidth=1.5) + contour!(axis, coordinate, coordinate, colorfield, + levels=[0.5], color=:white, linewidth=3) + scatter!(axis, particle_x, particle_y, color=(:black, 0.35), markersize=5) + arrows2d!(axis, interface_x, interface_y, normal_x, normal_y, + normalize=true, lengthscale=0.3, color=:black, + shaftwidth=3, tipwidth=14, tiplength=10) + xlims!(axis, extrema(coordinate)) + ylims!(axis, extrema(coordinate)) + Colorbar(fig[1, 2], heatmap, label="color field c") + + CairoMakie.save("colorfield_surface_normals.png", fig) +end +``` + +![Color-field level sets and interface-normal directions](colorfield_surface_normals.png) + +The particle phase has ``c \approx 1`` and the exterior has ``c \approx 0``. Consequently, +``\nabla c`` and the displayed unit normals point toward increasing ``c``. The arrows are +perpendicular to the color-field level sets; reversing the color convention reverses the arrows +without changing the interface geometry. #### Handling noise and errors in normal calculation -In regions distant from the interface, the calculated normals may be small or inaccurate due to the -smoothing kernel's support radius. To mitigate this: +Away from an interface, the exact color-field gradient vanishes, but particle disorder and +incomplete kernel support can produce small or poorly resolved normal estimates. To mitigate this: -1. Normals below a threshold are excluded from further calculations. -2. Curvature calculations use a corrected formulation to reduce errors near interface fringes. +1. Normals with insufficient particle support are discarded. +2. Morris-type calculations reject small gradients and can suppress normals in well-resolved + interior regions. +3. Curvature calculations use a corrected formulation to reduce errors near interface fringes. ```@autodocs Modules = [TrixiParticles] @@ -289,6 +386,17 @@ In the following table some values are shown for reference. The values marked wi | **Water** | 0.07288 [Lange](@cite Lange2005) | | **Mercury** | 0.486502 [Lange](@cite Lange2005) | +### Model configuration + +All surface tension coefficients must be finite and non-negative. A zero coefficient disables +the fluid-fluid surface force. Wall adhesion is controlled independently by the boundary's +`adhesion_coefficient`. + +`CohesionForceAkinci` only evaluates the pairwise cohesion and optional wall-adhesion forces. +It does not require surface normals or `reference_particle_spacing`. The full +`SurfaceTensionAkinci` model and both Morris models require a surface-normal method. When one +of these models is selected without an explicit method, `ColorfieldSurfaceNormal()` is used. + ### [Akinci-based intra-particle force surface tension and wall adhesion model](@id akinci_ipf) The [Akinci](@cite Akinci2013) model divides surface tension into distinct force components: @@ -350,6 +458,10 @@ A(r) = \frac{0.007}{h_c^{3.25}} \end{cases} ``` +The published adhesion kernel uses a three-dimensional normalization. In two-dimensional +simulations, `adhesion_coefficient` is therefore an empirical numerical parameter and may need +to be adjusted when changing the particle spacing or smoothing length. + --- ### [Morris surface tension model](@id morris_csf) 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/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..5556b13546 100644 --- a/src/schemes/fluid/entropically_damped_sph/system.jl +++ b/src/schemes/fluid/entropically_damped_sph/system.jl @@ -50,10 +50,13 @@ See [Entropically Damped Artificial Compressibility for SPH](@ref edac) for more The keyword argument `acceleration` should be used instead for gravity-like source terms. - `surface_tension`: Surface tension model used for this SPH system. (default: no surface tension) -- `surface_normal_method`: The surface normal method to be used for this SPH system. - (default: no surface normal method or `ColorfieldSurfaceNormal()` if a surface_tension model is used) +- `surface_normal_method`: Method used to estimate fluid-interface normals. This can be + configured independently for interface analysis and output and is + also used by models that require interface geometry. The default is + `nothing`; `ColorfieldSurfaceNormal()` is selected automatically + when the surface tension model requires normals. - `reference_particle_spacing`: The reference particle spacing used for weighting values at the boundary, - which currently is only needed when using surface tension. + which is needed when using a surface-normal method. - `color_value`: Integer label used for calculation of surface normals. Currently this is only used together with [`BoundaryModelDummyParticles`](@ref) and [`ColorfieldSurfaceNormal`](@ref): fluid-boundary normal evaluation @@ -119,12 +122,11 @@ function EntropicallyDampedSPHSystem(initial_condition; smoothing_kernel, smooth throw(ArgumentError("`acceleration` must be of length $NDIMS for a $(NDIMS)D problem")) end - if surface_tension !== nothing && surface_normal_method === nothing - surface_normal_method = ColorfieldSurfaceNormal() - end + surface_normal_method = default_surface_normal_method(surface_tension, + surface_normal_method) if surface_normal_method !== nothing && reference_particle_spacing < eps() - throw(ArgumentError("`reference_particle_spacing` must be set to a positive value when using `ColorfieldSurfaceNormal` or a surface tension model")) + throw(ArgumentError("`reference_particle_spacing` must be set to a positive value when using a surface-normal method")) end if correction isa ShepardKernelCorrection && @@ -308,7 +310,7 @@ function update_final!(system::EntropicallyDampedSPHSystem, v, u, v_ode, u_ode, kwargs...) (; surface_tension) = system - # Surface normal of neighbor and boundary needs to have been calculated already + # Surface-tension formulations using curvature or stress require previously computed normals. compute_curvature!(system, surface_tension, v, u, v_ode, u_ode, semi, t) compute_stress_tensors!(system, surface_tension, v, u, v_ode, u_ode, semi, t) update_average_pressure!(system, system.average_pressure_reduction, v_ode, u_ode, semi) 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..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. @@ -17,8 +19,22 @@ end function ColorfieldSurfaceNormal(; boundary_contact_threshold=0.1, interface_threshold=0.01, ideal_density_threshold=0.0) - return ColorfieldSurfaceNormal(boundary_contact_threshold, interface_threshold, - ideal_density_threshold) + thresholds = (boundary_contact_threshold, interface_threshold, + ideal_density_threshold) + if !all(threshold -> threshold isa Real && isfinite(threshold), thresholds) + throw(ArgumentError("surface-normal thresholds must be finite real numbers")) + end + + thresholds = promote(thresholds...) + return ColorfieldSurfaceNormal(thresholds...) +end + +@inline function default_surface_normal_method(surface_tension, surface_normal_method) + if isnothing(surface_normal_method) && requires_surface_normal(surface_tension) + return ColorfieldSurfaceNormal() + end + + return surface_normal_method end function create_cache_surface_normal(surface_normal_method, ELTYPE, NDIMS, nparticles) diff --git a/src/schemes/fluid/surface_tension.jl b/src/schemes/fluid/surface_tension.jl index 5656e95e12..fbb210e0c7 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 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 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 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 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) @@ -121,17 +149,24 @@ end pos_diff, distance) (; surface_tension_coefficient) = surface_tension - # Eq. 2 + distance >= support_radius && return zero(pos_diff) + + # 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 @@ -141,18 +176,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 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) + 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 # `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..4b0a1ff414 100644 --- a/src/schemes/fluid/weakly_compressible_sph/system.jl +++ b/src/schemes/fluid/weakly_compressible_sph/system.jl @@ -53,10 +53,13 @@ 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 a surface_tension model is used) +- `surface_normal_method`: Method used to estimate fluid-interface normals. This can be + configured independently for interface analysis and output and is + also used by models that require interface geometry. The default is + `nothing`; `ColorfieldSurfaceNormal()` is selected automatically + when the surface tension model requires normals. - `reference_particle_spacing`: The reference particle spacing used for weighting values at the boundary, - which currently is only needed when using surface tension. + which is needed when using a surface-normal method. - `color_value`: Integer label used for calculation of surface normals. Currently this is only used together with [`BoundaryModelDummyParticles`](@ref) and [`ColorfieldSurfaceNormal`](@ref): fluid-boundary normal evaluation @@ -130,12 +133,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, @@ -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 7fe8abbd97..687de166ef 100644 --- a/test/schemes/fluid/surface_tension.jl +++ b/test/schemes/fluid/surface_tension.jl @@ -1,5 +1,145 @@ - @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} + + 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 + 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, color_value=0) + edac = EntropicallyDampedSPHSystem(initial_condition; smoothing_kernel, + smoothing_length, sound_speed=10.0, + density_calculator=SummationDensity(), + surface_tension, color_value=0) + + 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) + + 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 isnothing(TrixiParticles.check_system_color((wcsph, edac))) + + @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, + 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 + 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 @@ -41,6 +181,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 @@ -88,6 +242,28 @@ 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) + + 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 "compute_stress_tensors! (MomentumMorris)" begin