diff --git a/NEWS.md b/NEWS.md index ad3e44c715..c294a476b1 100644 --- a/NEWS.md +++ b/NEWS.md @@ -4,9 +4,16 @@ TrixiParticles.jl follows the interpretation of [semantic versioning (semver)](https://julialang.github.io/Pkg.jl/dev/compatibility/#Version-specifier-format-1) used in the Julia ecosystem. Notable changes will be documented in this file for human readability. -## Version 0.5.3 - -### Features +## Version 0.6.0 + +### API Changes + +- Renamed the fluid-system keyword `surface_normal_method` to `surface_method` and added + detection-only surface methods. The old constructor keyword and accessor are deprecated. + +## Version 0.5.3 + +### Features - Added normal vectors to `InitialCondition`, with automatic computation for boundary particles generated by `RectangularTank` and `SphereShape` (#1036). diff --git a/docs/src/systems/fluid.md b/docs/src/systems/fluid.md index c9b9860ab9..45298f01a2 100644 --- a/docs/src/systems/fluid.md +++ b/docs/src/systems/fluid.md @@ -211,7 +211,7 @@ Pages = [joinpath("general", "corrections.jl")] --- -## [Surface Normals](@id surface_normal) +## [Surface Detection And Normals](@id surface_normal) ### Overview of surface normal calculation in SPH @@ -248,6 +248,28 @@ The calculated normals are normalized to unit vectors: Normalization ensures that the magnitude of the normals does not bias the curvature calculations or the resulting surface tension forces. +#### Surface methods and activity + +Fluid systems configure interface geometry with the `surface_method` keyword. Every surface +method computes a smooth `surface_activity`. Methods derived from +`AbstractSurfaceNormalMethod` additionally provide a normal, so normal calculation always +includes detection. + +`ColorfieldSurfaceDetection` computes activity only. `ColorfieldSurfaceNormal` uses the same +colorfield accumulation and additionally filters and stores the gradient as a surface normal. +Different `color_value`s detect interfaces between represented liquids. A constant nonzero +color detects a free surface because its kernel support ends at the unrepresented exterior. +Equal colors do not create an internal interface. + +`surface_activity` is available in particle VTK output and as a custom quantity for +`SolutionSavingCallback` and `PostprocessCallback`. `surf_normal` is written only for +normal-capable methods. + +Point and plane interpolation evaluate the same colorfield gradient. With `cut_off_bnd=true`, +kernel-weighted color contributions also determine whether a point belongs to the reference +phase. This prevents extrapolation through both free surfaces and interfaces with another +liquid while retaining the existing solid-boundary cutoff. + #### Handling noise and errors in normal calculation In regions distant from the interface, the calculated normals may be small or inaccurate due to the diff --git a/examples/fluid/sphere_surface_tension_2d.jl b/examples/fluid/sphere_surface_tension_2d.jl index 1abcdcb149..b22e519e1d 100644 --- a/examples/fluid/sphere_surface_tension_2d.jl +++ b/examples/fluid/sphere_surface_tension_2d.jl @@ -53,7 +53,7 @@ fluid_system = EntropicallyDampedSPHSystem(fluid; smoothing_kernel=fluid_smoothi density_calculator=ContinuityDensity(), reference_particle_spacing=particle_spacing, acceleration=zeros(length(fluid_size)), - surface_normal_method=ColorfieldSurfaceNormal(), + surface_method=ColorfieldSurfaceNormal(), surface_tension=SurfaceTensionMorris(surface_tension_coefficient=50 * 0.0728)) diff --git a/src/TrixiParticles.jl b/src/TrixiParticles.jl index ee3bf532b5..3b9f23d91d 100644 --- a/src/TrixiParticles.jl +++ b/src/TrixiParticles.jl @@ -115,7 +115,9 @@ export interpolate_line, interpolate_points, interpolate_plane_3d, interpolate_p interpolate_plane_2d_vtk export SurfaceTensionAkinci, CohesionForceAkinci, SurfaceTensionMorris, SurfaceTensionMomentumMorris -export ColorfieldSurfaceNormal +export AbstractSurfaceMethod, AbstractSurfaceNormalMethod, ColorfieldSurfaceDetection, + ColorfieldSurfaceNormal, surface_method, computes_surface_normal, surface_activity, + surface_normal export SymplecticPositionVerlet export coordinates_eltype diff --git a/src/general/custom_quantities.jl b/src/general/custom_quantities.jl index 3f1f74c637..0715678c2a 100644 --- a/src/general/custom_quantities.jl +++ b/src/general/custom_quantities.jl @@ -141,3 +141,24 @@ end function avg_density(system, dv_ode, du_ode, v_ode, u_ode, semi, t) return NaN end + +""" + surface_activity + +Return the per-particle smooth surface activity, or `nothing` for systems without a +configured surface method. +""" +function surface_activity(system::AbstractFluidSystem, dv_ode, du_ode, v_ode, u_ode, + semi, t) + hasproperty(system.cache, :surface_activity) || return nothing + return view(system.cache.surface_activity, each_active_particle(system)) +end + +surface_activity(system, dv_ode, du_ode, v_ode, u_ode, semi, t) = nothing + +function surface_normal(system::AbstractFluidSystem, dv_ode, du_ode, v_ode, u_ode, semi, t) + computes_surface_normal(surface_method(system)) || return nothing + return view(system.cache.surface_normal, :, each_active_particle(system)) +end + +surface_normal(system, dv_ode, du_ode, v_ode, u_ode, semi, t) = nothing diff --git a/src/general/interpolation.jl b/src/general/interpolation.jl index 961224fab7..75c5f9e777 100644 --- a/src/general/interpolation.jl +++ b/src/general/interpolation.jl @@ -23,11 +23,11 @@ See also: [`interpolate_plane_2d_vtk`](@ref), [`interpolate_plane_3d`](@ref), # Keywords - `smoothing_length=initial_smoothing_length(ref_system)`: The smoothing length used in the interpolation. -- `cut_off_bnd=true`: Boolean to indicate if quantities should be set to `NaN` when the point - is "closer" to the boundary than to the fluid in a kernel-weighted sense. - Or, in more detail, when the boundary has more influence than the fluid - on the density summation in this point, i.e., when the boundary particles - add more kernel-weighted mass than the fluid particles. +- `cut_off_bnd=true`: Boolean to indicate if quantities should be set to `NaN` outside the + interpolated fluid domain. Boundaries are detected by comparing their + kernel-weighted mass to that of the fluid. When `ref_system` uses + a colorfield surface method, free and multiphase surfaces are detected + from the same color contributions used for particle surface detection. - `clip_negative_pressure=false`: One common approach in SPH models is to clip negative pressure values, but this is unphysical. Instead we clip here during interpolation thus only impacting the local interpolated value. @@ -41,7 +41,7 @@ See also: [`interpolate_plane_2d_vtk`](@ref), [`interpolate_plane_3d`](@ref), !!! note - The interpolation accuracy is subject to the density of particles and the chosen smoothing length. - - With `cut_off_bnd`, a density-based estimation of the surface is used, which is not as + - With `cut_off_bnd`, a kernel-based estimation of the surface is used, which is not as accurate as a real surface reconstruction. # Examples @@ -115,11 +115,9 @@ See also: [`interpolate_plane_2d`](@ref), [`interpolate_plane_3d`](@ref), - `smoothing_length=initial_smoothing_length(ref_system)`: The smoothing length used in the interpolation. - `output_directory="out"`: Directory to save the VTI file. - `filename="plane"`: Name of the VTI file. -- `cut_off_bnd=true`: Boolean to indicate if quantities should be set to `NaN` when the point - is "closer" to the boundary than to the fluid in a kernel-weighted sense. - Or, in more detail, when the boundary has more influence than the fluid - on the density summation in this point, i.e., when the boundary particles - add more kernel-weighted mass than the fluid particles. +- `cut_off_bnd=true`: Set quantities outside the interpolated fluid domain to `NaN`. + Solid boundaries use kernel-weighted mass. A configured colorfield + surface method additionally uses phase-color contributions. - `clip_negative_pressure=false`: One common approach in SPH models is to clip negative pressure values, but this is unphysical. Instead we clip here during interpolation thus only impacting the local interpolated value. @@ -130,7 +128,7 @@ See also: [`interpolate_plane_2d`](@ref), [`interpolate_plane_3d`](@ref), !!! note - The interpolation accuracy is subject to the density of particles and the chosen smoothing length. - - With `cut_off_bnd`, a density-based estimation of the surface is used, which is not as + - With `cut_off_bnd`, a kernel-based estimation of the surface is used, which is not as accurate as a real surface reconstruction. # Examples @@ -177,6 +175,10 @@ function interpolate_plane_2d_vtk(min_corner, max_corner, resolution, semi, ref_ vtk["density"] = density vtk["velocity"] = velocity vtk["pressure"] = pressure + if hasproperty(results, :surface_activity) + vtk["surface_activity"] = reshape(results.surface_activity, + length(x_range), length(y_range)) + end end end @@ -250,11 +252,9 @@ See also: [`interpolate_plane_2d`](@ref), [`interpolate_plane_2d_vtk`](@ref), # Keywords - `smoothing_length=initial_smoothing_length(ref_system)`: The smoothing length used in the interpolation. -- `cut_off_bnd=true`: Boolean to indicate if quantities should be set to `NaN` when the point - is "closer" to the boundary than to the fluid in a kernel-weighted sense. - Or, in more detail, when the boundary has more influence than the fluid - on the density summation in this point, i.e., when the boundary particles - add more kernel-weighted mass than the fluid particles. +- `cut_off_bnd=true`: Set quantities outside the interpolated fluid domain to `NaN`. + Solid boundaries use kernel-weighted mass. A configured colorfield + surface method additionally uses phase-color contributions. - `clip_negative_pressure=false`: One common approach in SPH models is to clip negative pressure values, but this is unphysical. Instead we clip here during interpolation thus only impacting the local interpolated value. @@ -268,7 +268,7 @@ See also: [`interpolate_plane_2d`](@ref), [`interpolate_plane_2d_vtk`](@ref), !!! note - The interpolation accuracy is subject to the density of particles and the chosen smoothing length. - - With `cut_off_bnd`, a density-based estimation of the surface is used which is not as + - With `cut_off_bnd`, a kernel-based estimation of the surface is used which is not as accurate as a real surface reconstruction. # Examples @@ -346,11 +346,9 @@ See also: [`interpolate_points`](@ref), [`interpolate_plane_2d`](@ref), # Keywords - `endpoint=true`: A boolean to include (`true`) or exclude (`false`) the end point in the interpolation. - `smoothing_length=initial_smoothing_length(ref_system)`: The smoothing length used in the interpolation. -- `cut_off_bnd=true`: Boolean to indicate if quantities should be set to `NaN` when the point - is "closer" to the boundary than to the fluid in a kernel-weighted sense. - Or, in more detail, when the boundary has more influence than the fluid - on the density summation in this point, i.e., when the boundary particles - add more kernel-weighted mass than the fluid particles. +- `cut_off_bnd=true`: Set quantities outside the interpolated fluid domain to `NaN`. + Solid boundaries use kernel-weighted mass. A configured colorfield + surface method additionally uses phase-color contributions. - `clip_negative_pressure=false`: One common approach in SPH models is to clip negative pressure values, but this is unphysical. Instead we clip here during interpolation thus only impacting the local interpolated value. @@ -366,7 +364,7 @@ See also: [`interpolate_points`](@ref), [`interpolate_plane_2d`](@ref), - This function is particularly useful for analyzing gradients or creating visualizations along a specified line in the SPH simulation domain. - The interpolation accuracy is subject to the density of particles and the chosen smoothing length. - - With `cut_off_bnd`, a density-based estimation of the surface is used which is not as + - With `cut_off_bnd`, a kernel-based estimation of the surface is used which is not as accurate as a real surface reconstruction. # Examples @@ -431,11 +429,9 @@ See also: [`interpolate_line`](@ref), [`interpolate_plane_2d`](@ref), # Keywords - `smoothing_length=initial_smoothing_length(ref_system)`: The smoothing length used in the interpolation. -- `cut_off_bnd=true`: Boolean to indicate if quantities should be set to `NaN` when the point - is "closer" to the boundary than to the fluid in a kernel-weighted sense. - Or, in more detail, when the boundary has more influence than the fluid - on the density summation in this point, i.e., when the boundary particles - add more kernel-weighted mass than the fluid particles. +- `cut_off_bnd=true`: Set quantities outside the interpolated fluid domain to `NaN`. + Solid boundaries use kernel-weighted mass. A configured colorfield + surface method additionally uses phase-color contributions. - `clip_negative_pressure=false`: One common approach in SPH models is to clip negative pressure values, but this is unphysical. Instead we clip here during interpolation thus only impacting the local interpolated value. @@ -463,7 +459,7 @@ results = interpolate_points(points, semi, ref_system, sol) - This function is particularly useful for analyzing gradients or creating visualizations along a specified line in the SPH simulation domain. - The interpolation accuracy is subject to the density of particles and the chosen smoothing length. - - With `cut_off_bnd`, a density-based estimation of the surface is used which is not as + - With `cut_off_bnd`, a kernel-based estimation of the surface is used which is not as accurate as a real surface reconstruction. """ @inline function interpolate_points(point_coords, semi, ref_system, sol::ODESolution; @@ -555,10 +551,22 @@ end n_points = size(point_coords, 2) ELTYPE = eltype(point_coords) + surface_method_ = surface_method(ref_system) + detect_surface = is_colorfield_surface_method(surface_method_) computed_density = allocate(semi.parallelization_backend, ELTYPE, n_points) other_density = allocate(semi.parallelization_backend, ELTYPE, n_points) shepard_coefficient = allocate(semi.parallelization_backend, ELTYPE, n_points) neighbor_count = allocate(semi.parallelization_backend, Int, n_points) + surface_gradient = detect_surface ? + allocate(semi.parallelization_backend, ELTYPE, + (ndims(ref_system), n_points)) : nothing + surface_activity_ = detect_surface ? + allocate(semi.parallelization_backend, ELTYPE, n_points) : nothing + reference_colorfield = detect_surface ? + allocate(semi.parallelization_backend, ELTYPE, n_points) : + nothing + other_colorfield = detect_surface ? + allocate(semi.parallelization_backend, ELTYPE, n_points) : nothing # The wall velocity considers more neighbors, so we need to use # a different Shepard coefficient. shepard_coefficient_wall = allocate(semi.parallelization_backend, ELTYPE, n_points) @@ -568,15 +576,26 @@ end set_zero!(other_density) set_zero!(shepard_coefficient) set_zero!(neighbor_count) + if detect_surface + set_zero!(surface_gradient) + set_zero!(surface_activity_) + set_zero!(reference_colorfield) + set_zero!(other_colorfield) + end cache = create_cache_interpolation(ref_system, n_points, semi) ref_id = system_indices(ref_system, semi) ref_smoothing_kernel = ref_system.smoothing_kernel + interpolation_surface_threshold = detect_surface ? + surface_method_.interpolation_surface_threshold : + zero(ELTYPE) + reference_color = detect_surface ? ref_system.cache.color : 0 # If we neither cut off at the boundary nor include the boundary wall velocity, # we only need to iterate over the reference system. - systems = (cut_off_bnd || include_wall_velocity) ? semi : (ref_system,) + systems = (cut_off_bnd || include_wall_velocity || detect_surface) ? + semi : (ref_system,) foreach_system(systems) do neighbor_system system_id = system_indices(neighbor_system, semi) @@ -586,6 +605,12 @@ end u = wrap_u(u_ode, neighbor_system, semi) neighbor_coords = current_coordinates(u, neighbor_system) + contributes_surface = detect_surface && + has_system_interaction(ref_system, neighbor_system, semi) && + (contributes_to_colorfield(neighbor_system) || + contributes_boundary_colorfield(neighbor_system)) + surface_color = contributes_to_colorfield(neighbor_system) ? + neighbor_system.cache.color : reference_color foreach_point_neighbor(point_coords, neighbor_coords, nhs; parallelization_backend) do point, neighbor, pos_diff, @@ -594,6 +619,22 @@ end volume_b = m_b / current_density(v, neighbor_system, neighbor) W_ab = kernel(ref_smoothing_kernel, distance, smoothing_length) + if contributes_surface + grad_kernel = kernel_grad(ref_smoothing_kernel, pos_diff, distance, + smoothing_length) + for i in 1:ndims(ref_system) + surface_gradient[i, point] += volume_b * surface_color * grad_kernel[i] + end + + if neighbor_system isa AbstractFluidSystem + if surface_color == reference_color + reference_colorfield[point] += volume_b * W_ab + else + other_colorfield[point] += volume_b * W_ab + end + end + end + if include_wall_velocity # The wall velocity considers more neighbors, so we need to use # a different Shepard coefficient. @@ -610,7 +651,10 @@ end interpolate_system!(cache, v, neighbor_system, point, neighbor, volume_b, W_ab, clip_negative_pressure) else - other_density[point] += m_b * W_ab + if cut_off_bnd && + (!detect_surface || !(neighbor_system isa AbstractFluidSystem)) + other_density[point] += m_b * W_ab + end if include_wall_velocity velocity_neighbor_ = current_velocity(v, neighbor_system, neighbor) @@ -622,17 +666,38 @@ end end end - neighbor_count[point] += 1 + if system_id == ref_id || cut_off_bnd || include_wall_velocity + neighbor_count[point] += 1 + end end end @threaded parallelization_backend for point in axes(point_coords, 2) + if detect_surface + normal_norm = zero(ELTYPE) + for i in 1:ndims(ref_system) + normal_norm += surface_gradient[i, point]^2 + end + surface_activity_[point] = gradient_surface_activity(sqrt(normal_norm), + compact_support(ref_smoothing_kernel, + smoothing_length), + surface_method_) + end + + outside_reference_phase = detect_surface && + (reference_colorfield[point] < + interpolation_surface_threshold || + other_colorfield[point] > reference_colorfield[point]) cut_off = computed_density[point] < eps() || - (cut_off_bnd && other_density[point] > computed_density[point]) + (cut_off_bnd && (other_density[point] > computed_density[point] || + outside_reference_phase)) if cut_off # Return NaN values that can be filtered out in ParaView computed_density[point] = NaN neighbor_count[point] = 0 + if detect_surface + surface_activity_[point] = NaN + end # We need to convert the `NamedTuple` to a `Tuple` for GPU compatibility foreach(Tuple(cache)) do field @@ -658,7 +723,10 @@ end end end - return (; computed_density, point_coords, neighbor_count, cache...) + surface_detection_output = detect_surface ? + (; surface_activity=surface_activity_) : (;) + return (; computed_density, point_coords, neighbor_count, surface_detection_output..., + cache...) end @inline function create_cache_interpolation(ref_system::AbstractFluidSystem, n_points, semi) diff --git a/src/general/semidiscretization.jl b/src/general/semidiscretization.jl index a9ed769dc6..4ef14ebf7c 100644 --- a/src/general/semidiscretization.jl +++ b/src/general/semidiscretization.jl @@ -920,22 +920,15 @@ 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 is_colorfield_surface_method(surface_method(system)) end if requires_color_check + system_ids = findall(system -> system isa AbstractFluidSystem && + !(system isa ParticlePackingSystem), systems) - # Systems that contribute to the colorfield/contact logic. - system_ids = findall(system -> (system isa AbstractFluidSystem && - !(system isa ParticlePackingSystem)) || - system isa WallBoundarySystem || - system isa - RigidBodySystem{<:BoundaryModelDummyParticles}, - systems) - - if length(system_ids) > 1 && sum(i -> systems[i].cache.color, system_ids) == 0 - throw(ArgumentError("If `ColorfieldSurfaceNormal` or a surface tension model is used, at least one participating system must have a color different from 0.")) + if all(i -> iszero(systems[i].cache.color), system_ids) + throw(ArgumentError("If a colorfield surface method or a surface tension model is used, at least one participating system must have a color different from 0.")) end end end diff --git a/src/io/io.jl b/src/io/io.jl index 692dfd8528..d1b3e142fa 100644 --- a/src/io/io.jl +++ b/src/io/io.jl @@ -86,7 +86,7 @@ function add_system_data!(system_data, system::AbstractFluidSystem) system_data["pressure_acceleration_formulation"] = nameof(system.pressure_acceleration_formulation) add_system_data!(system_data, shifting_technique(system)) add_system_data!(system_data, system.surface_tension) - add_system_data!(system_data, system.surface_normal_method) + add_system_data!(system_data, system.surface_method) add_system_data!(system_data, system.viscosity) add_system_data!(system_data, system.correction) add_system_data!(system_data, system_state_equation(system)) @@ -107,6 +107,7 @@ function add_system_data!(system_data, system::ImplicitIncompressibleSPHSystem) system_data["acceleration"] = system.acceleration system_data["pressure_acceleration_formulation"] = nameof(system.pressure_acceleration_formulation) add_system_data!(system_data, shifting_technique(system)) + add_system_data!(system_data, system.surface_method) add_system_data!(system_data, system.viscosity) end @@ -315,11 +316,15 @@ function add_system_data!(system_data, system_data["surface_tension"]["surface_tension_coefficient"] = surface_tension.surface_tension_coefficient end -function add_system_data!(system_data, surface_normal_method::ColorfieldSurfaceNormal) - system_data["surface_normal_method"] = Dict{String, Any}() - system_data["surface_normal_method"]["model"] = type2string(surface_normal_method) - system_data["surface_normal_method"]["boundary_contact_threshold"] = surface_normal_method.boundary_contact_threshold - system_data["surface_normal_method"]["ideal_density_threshold"] = surface_normal_method.ideal_density_threshold +function add_system_data!(system_data, surface_method_::ColorfieldSurfaceMethod) + system_data["surface_method"] = Dict{String, Any}() + system_data["surface_method"]["model"] = type2string(surface_method_) + system_data["surface_method"]["computes_surface_normal"] = computes_surface_normal(surface_method_) + system_data["surface_method"]["boundary_contact_threshold"] = surface_method_.boundary_contact_threshold + system_data["surface_method"]["interface_threshold"] = surface_method_.interface_threshold + system_data["surface_method"]["ideal_density_threshold"] = surface_method_.ideal_density_threshold + system_data["surface_method"]["interface_taper_start"] = surface_method_.interface_taper_start + system_data["surface_method"]["interpolation_surface_threshold"] = surface_method_.interpolation_surface_threshold end function add_system_data!(system_data, boundary_zone::BoundaryZone, indice) diff --git a/src/io/write_vtk.jl b/src/io/write_vtk.jl index f242fd74b2..43320b2e4e 100644 --- a/src/io/write_vtk.jl +++ b/src/io/write_vtk.jl @@ -329,11 +329,15 @@ function write2vtk!(vtk, v, u, t, system::AbstractFluidSystem) vtk["pressure"] = [current_pressure(v, system, particle) for particle in eachparticle(system)] - if system.surface_normal_method !== nothing - vtk["surf_normal"] = [surface_normal(system, particle) - for particle in eachparticle(system)] + if system.surface_method !== nothing + vtk["surface_activity"] = system.cache.surface_activity vtk["neighbor_count"] = system.cache.neighbor_count vtk["color"] = system.cache.color + + if computes_surface_normal(system.surface_method) + vtk["surf_normal"] = [surface_normal(system, particle) + for particle in eachparticle(system)] + end end if system.surface_tension isa SurfaceTensionMorris || diff --git a/src/schemes/boundary/wall_boundary/system.jl b/src/schemes/boundary/wall_boundary/system.jl index f9864ecc05..3791970de4 100644 --- a/src/schemes/boundary/wall_boundary/system.jl +++ b/src/schemes/boundary/wall_boundary/system.jl @@ -14,9 +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. +- `color_value`: Integer label used for colorfield surface calculations. Currently this is only used together with [`BoundaryModelDummyParticles`](@ref) and - [`ColorfieldSurfaceNormal`](@ref): fluid-boundary normal evaluation + colorfield surface methods: fluid-boundary surface evaluation reads the resulting boundary colorfield to detect wall contact. """ struct WallBoundarySystem{BM, ELTYPE <: Real, NDIMS, IC, CO, M, IM, diff --git a/src/schemes/fluid/entropically_damped_sph/system.jl b/src/schemes/fluid/entropically_damped_sph/system.jl index 6e03cf7837..e8386388d9 100644 --- a/src/schemes/fluid/entropically_damped_sph/system.jl +++ b/src/schemes/fluid/entropically_damped_sph/system.jl @@ -6,7 +6,7 @@ shifting_technique=nothing, alpha=0.5, viscosity=nothing, acceleration=ntuple(_ -> 0.0, NDIMS), surface_tension=nothing, - surface_normal_method=nothing, buffer_size=nothing, + surface_method=nothing, buffer_size=nothing, reference_particle_spacing=0.0, color_value=1, source_terms=nothing) @@ -50,14 +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) -- `reference_particle_spacing`: The reference particle spacing used for weighting values at the boundary, - which currently is only needed when using surface tension. -- `color_value`: Integer label used for calculation of surface normals. - Currently this is only used together with [`BoundaryModelDummyParticles`](@ref) and - [`ColorfieldSurfaceNormal`](@ref): fluid-boundary normal evaluation - reads the resulting boundary colorfield to detect wall contact. +- `surface_method`: Surface detection or normal method used by this system. + Methods that compute normals always also compute surface + activity. The default is `nothing`, or + `ColorfieldSurfaceNormal()` when required by surface tension. +- `reference_particle_spacing`: Reference spacing required by colorfield surface methods. +- `color_value`: Scalar contributed to colorfield surface detection. Different + values identify represented fluid-fluid interfaces. """ struct EntropicallyDampedSPHSystem{NDIMS, ELTYPE <: Real, IC, M, DC, K, V, COR, PF, TV, @@ -77,7 +76,7 @@ struct EntropicallyDampedSPHSystem{NDIMS, ELTYPE <: Real, IC, M, DC, K, V, COR, average_pressure_reduction :: AVGP source_terms :: ST surface_tension :: SRFT - surface_normal_method :: SRFN + surface_method :: SRFN buffer :: B particle_refinement :: PR cache :: C @@ -95,7 +94,8 @@ function EntropicallyDampedSPHSystem(initial_condition; smoothing_kernel, smooth ndims(smoothing_kernel)), correction=nothing, source_terms=nothing, surface_tension=nothing, - surface_normal_method=nothing, buffer_size=nothing, + surface_method=nothing, surface_normal_method=nothing, + buffer_size=nothing, reference_particle_spacing=0.0, color_value=1) buffer = isnothing(buffer_size) ? nothing : SystemBuffer(nparticles(initial_condition), buffer_size) @@ -119,12 +119,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_method = select_surface_method(surface_tension, surface_method, + 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")) + if is_colorfield_surface_method(surface_method) && reference_particle_spacing < eps() + throw(ArgumentError("`reference_particle_spacing` must be set to a positive value when using a colorfield surface method")) end if correction isa ShepardKernelCorrection && @@ -145,8 +144,7 @@ function EntropicallyDampedSPHSystem(initial_condition; smoothing_kernel, smooth create_cache_shifting(initial_condition, shifting_technique)..., create_cache_avg_pressure_reduction(initial_condition, avg_pressure_reduction)..., - create_cache_surface_normal(surface_normal_method, ELTYPE, NDIMS, - n_particles)..., + create_cache_surface(surface_method, ELTYPE, NDIMS, n_particles)..., create_cache_surface_tension(surface_tension, ELTYPE, NDIMS, n_particles)..., create_cache_refinement(initial_condition, particle_refinement, @@ -168,7 +166,7 @@ function EntropicallyDampedSPHSystem(initial_condition; smoothing_kernel, smooth typeof(viscosity), typeof(correction), typeof(pressure_acceleration), typeof(shifting_technique), typeof(avg_pressure_reduction), typeof(source_terms), - typeof(surface_tension), typeof(surface_normal_method), + typeof(surface_tension), typeof(surface_method), typeof(buffer), Nothing, typeof(cache)}(initial_condition, mass, density_calculator, smoothing_kernel, sound_speed, viscosity, @@ -176,7 +174,7 @@ function EntropicallyDampedSPHSystem(initial_condition; smoothing_kernel, smooth pressure_acceleration, shifting_technique, avg_pressure_reduction, source_terms, surface_tension, - surface_normal_method, buffer, + surface_method, buffer, particle_refinement, cache) end @@ -199,7 +197,7 @@ function Base.show(io::IO, system::EntropicallyDampedSPHSystem) print(io, ", ", system.smoothing_kernel) print(io, ", ", system.acceleration) print(io, ", ", system.surface_tension) - print(io, ", ", system.surface_normal_method) + print(io, ", ", system.surface_method) print(io, ") with ", nparticles(system), " particles") end @@ -228,7 +226,7 @@ function Base.show(io::IO, ::MIME"text/plain", system::EntropicallyDampedSPHSyst typeof(system.average_pressure_reduction).parameters[1] ? "yes" : "no") summary_line(io, "acceleration", system.acceleration) summary_line(io, "surface tension", system.surface_tension) - summary_line(io, "surface normal method", system.surface_normal_method) + summary_line(io, "surface method", system.surface_method) summary_footer(io) end end @@ -299,8 +297,7 @@ function update_quantities!(system::EntropicallyDampedSPHSystem, v, u, end function update_pressure!(system::EntropicallyDampedSPHSystem, v, u, v_ode, u_ode, semi, t) - compute_surface_normal!(system, system.surface_normal_method, v, u, v_ode, u_ode, semi, - t) + compute_surface!(system, system.surface_method, v, u, v_ode, u_ode, semi, t) compute_surface_delta_function!(system, system.surface_tension, semi) end diff --git a/src/schemes/fluid/fluid.jl b/src/schemes/fluid/fluid.jl index fb49a65359..317aff3b9c 100644 --- a/src/schemes/fluid/fluid.jl +++ b/src/schemes/fluid/fluid.jl @@ -246,14 +246,27 @@ end return nothing end -@inline function surface_normal_method(system::AbstractFluidSystem) - return system.surface_normal_method +@inline function surface_method(system::AbstractFluidSystem) + hasproperty(system, :surface_method) || return nothing + return system.surface_method end -@inline function surface_normal_method(system) +@inline function surface_method(system) return nothing end +function surface_normal_method(system) + Base.depwarn("`surface_normal_method(system)` is deprecated; use `surface_method(system)`", + :surface_normal_method) + method = surface_method(system) + return computes_surface_normal(method) ? method : nothing +end + +@inline contributes_to_colorfield(system) = false +@inline function contributes_to_colorfield(system::AbstractFluidSystem) + return hasproperty(system.cache, :color) +end + function restart_u(system::AbstractFluidSystem, data) inactive_coords = convert(coordinates_eltype(system), 1e16) coords_total = fill(inactive_coords, u_nvariables(system), @@ -297,13 +310,13 @@ 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 requires_surface_normal(fluid_system.surface_tension) foreach_system(systems) do neighbor - if neighbor isa AbstractFluidSystem && - isnothing(fluid_system.surface_tension) && - isnothing(fluid_system.surface_normal_method) - throw(ArgumentError("either none or all fluid systems in a simulation need " * - "to use a surface tension model or a surface normal method.")) + if neighbor isa AbstractFluidSystem && !(neighbor isa ParticlePackingSystem) && + !computes_surface_normal(surface_method(neighbor)) + throw(ArgumentError("all interacting fluid systems must use a surface method " * + "that computes normals when a surface-tension model " * + "requires interface normals")) end end end diff --git a/src/schemes/fluid/implicit_incompressible_sph/system.jl b/src/schemes/fluid/implicit_incompressible_sph/system.jl index e4ef03584f..065f74bf29 100644 --- a/src/schemes/fluid/implicit_incompressible_sph/system.jl +++ b/src/schemes/fluid/implicit_incompressible_sph/system.jl @@ -4,7 +4,9 @@ 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, + surface_method=nothing, + reference_particle_spacing=0.0, color_value=1) System for particles of a fluid. The system employs implicit incompressible SPH (IISPH), iteratively solving a linear system @@ -30,9 +32,12 @@ 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 +- `surface_method`: Optional surface detection or normal method. +- `reference_particle_spacing`: Reference spacing required by colorfield surface methods. +- `color_value`: Scalar contributed to colorfield surface calculations. """ struct ImplicitIncompressibleSPHSystem{NDIMS, ELTYPE <: Real, ARRAY1D, ARRAY2D, - IC, K, V, PF, C} <: AbstractFluidSystem{NDIMS} + IC, K, V, PF, SM, C} <: AbstractFluidSystem{NDIMS} initial_condition :: IC mass :: ARRAY1D # Array{ELTYPE, 1} pressure :: ARRAY1D @@ -42,7 +47,7 @@ struct ImplicitIncompressibleSPHSystem{NDIMS, ELTYPE <: Real, ARRAY1D, ARRAY2D, acceleration :: SVector{NDIMS, ELTYPE} viscosity :: V pressure_acceleration_formulation :: PF - surface_normal_method :: Nothing # TODO + surface_method :: SM surface_tension :: Nothing # TODO particle_refinement :: Nothing # TODO density :: ARRAY1D @@ -71,7 +76,10 @@ 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, + surface_method=nothing, + surface_normal_method=nothing, + reference_particle_spacing=0.0, color_value=1) particle_refinement = nothing # TODO surface_tension = nothing # TODO @@ -112,6 +120,12 @@ function ImplicitIncompressibleSPHSystem(initial_condition; smoothing_kernel, throw(ArgumentError("`time_step` must be a positive number")) end + surface_method = select_surface_method(surface_tension, surface_method, + surface_normal_method) + if is_colorfield_surface_method(surface_method) && reference_particle_spacing < eps() + throw(ArgumentError("`reference_particle_spacing` must be set to a positive value when using a colorfield surface method")) + end + pressure_acceleration = pressure_acceleration_summation_density density = copy(initial_condition.density) @@ -124,13 +138,19 @@ function ImplicitIncompressibleSPHSystem(initial_condition; smoothing_kernel, density_error = zeros(ELTYPE, n_particles) cache = (; + create_cache_surface(surface_method, ELTYPE, NDIMS, n_particles)..., create_cache_refinement(initial_condition, particle_refinement, - smoothing_length)...,) + smoothing_length)..., + color=Int(color_value)) + if reference_particle_spacing > 0 + cache = (; cache..., reference_particle_spacing) + end return ImplicitIncompressibleSPHSystem(initial_condition, mass, pressure, smoothing_kernel, smoothing_length, reference_density, acceleration_, viscosity, - pressure_acceleration, nothing, surface_tension, + pressure_acceleration, surface_method, + surface_tension, particle_refinement, density, predicted_density, advection_velocity, d_ii, a_ii, sum_d_ij_pj, sum_term, density_error, omega, max_error, @@ -216,6 +236,12 @@ function update_quantities!(system::ImplicitIncompressibleSPHSystem, v, u, semi) end +function update_pressure!(system::ImplicitIncompressibleSPHSystem, v, u, v_ode, u_ode, + semi, t) + compute_surface!(system, system.surface_method, v, u, v_ode, u_ode, semi, t) + return system +end + function update_implicit_sph!(semi, v_ode, u_ode, t) # This check is performed statically by the compiler and has no overhead if !any(system -> system isa ImplicitIncompressibleSPHSystem, semi.systems) diff --git a/src/schemes/fluid/surface_normal_sph.jl b/src/schemes/fluid/surface_normal_sph.jl index adbc9d7dbe..5d394c762f 100644 --- a/src/schemes/fluid/surface_normal_sph.jl +++ b/src/schemes/fluid/surface_normal_sph.jl @@ -1,58 +1,235 @@ +abstract type AbstractSurfaceMethod end +abstract type AbstractSurfaceNormalMethod <: AbstractSurfaceMethod end + @doc raw""" - ColorfieldSurfaceNormal(; boundary_contact_threshold=0.1, interface_threshold=0.01, - ideal_density_threshold=0.0) + ColorfieldSurfaceDetection(; boundary_contact_threshold=0.1, + interface_threshold=0.01, + ideal_density_threshold=0.0, + interface_taper_start=0.8, + interpolation_surface_threshold=0.45) + +Detect fluid surfaces from the magnitude of a colorfield gradient. Different +`color_value`s detect represented fluid-fluid interfaces, while incomplete constant-color +support detects a free surface. This method computes [`surface_activity`](@ref), but does +not expose a surface normal. +""" +struct ColorfieldSurfaceDetection{ELTYPE} <: AbstractSurfaceMethod + boundary_contact_threshold::ELTYPE + interface_threshold::ELTYPE + ideal_density_threshold::ELTYPE + interface_taper_start::ELTYPE + interpolation_surface_threshold::ELTYPE +end -Color field based computation of the interface normals. +@doc raw""" + ColorfieldSurfaceNormal(; boundary_contact_threshold=0.1, + interface_threshold=0.01, + ideal_density_threshold=0.0, + interface_taper_start=0.8, + interpolation_surface_threshold=0.45) + +Compute colorfield surface normals and [`surface_activity`](@ref). The detection stage is +identical to [`ColorfieldSurfaceDetection`](@ref). The raw gradient is filtered after its +magnitude has been converted to activity and is normalized when required by the configured +surface-tension model. # 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]` used to detect contact with + dummy-particle boundaries. +- `interface_threshold=0.01`: Finite, non-negative dimensionless gradient threshold. +- `ideal_density_threshold=0.0`: Optional neighbor-count heuristic for unrepresented exterior + phases. Zero disables it; keep it disabled for fully supported multiphase interfaces. +- `interface_taper_start=0.8`: Start of the smooth activity transition as a fraction of + `interface_threshold`. +- `interpolation_surface_threshold=0.45`: Minimum normalized reference-color contribution + retained by interpolated output. """ -struct ColorfieldSurfaceNormal{ELTYPE} +struct ColorfieldSurfaceNormal{ELTYPE} <: AbstractSurfaceNormalMethod boundary_contact_threshold::ELTYPE interface_threshold::ELTYPE ideal_density_threshold::ELTYPE + interface_taper_start::ELTYPE + interpolation_surface_threshold::ELTYPE end -function ColorfieldSurfaceNormal(; boundary_contact_threshold=0.1, interface_threshold=0.01, - ideal_density_threshold=0.0) - return ColorfieldSurfaceNormal(boundary_contact_threshold, interface_threshold, +const ColorfieldSurfaceMethod = Union{ColorfieldSurfaceDetection, ColorfieldSurfaceNormal} + +function validate_surface_threshold(threshold, name; upper_bound=nothing, + strict_upper_bound=false) + interval = if isnothing(upper_bound) + "non-negative" + elseif strict_upper_bound + "in [0, $upper_bound)" + else + "in [0, $upper_bound]" + end + threshold isa Real || + throw(ArgumentError("`$name` must be a finite real number $interval")) + + valid_upper_bound = isnothing(upper_bound) || + (strict_upper_bound ? threshold < upper_bound : + threshold <= upper_bound) + if !isfinite(threshold) || threshold < 0 || !valid_upper_bound + throw(ArgumentError("`$name` must be a finite real number $interval")) + end + + return threshold +end + +function colorfield_surface_parameters(; boundary_contact_threshold=0.1, + interface_threshold=0.01, + ideal_density_threshold=0.0, + interface_taper_start=0.8, + interpolation_surface_threshold=0.45) + boundary_threshold = validate_surface_threshold(boundary_contact_threshold, + "boundary_contact_threshold"; + upper_bound=1) + normal_threshold = validate_surface_threshold(interface_threshold, + "interface_threshold") + density_threshold = validate_surface_threshold(ideal_density_threshold, + "ideal_density_threshold"; + upper_bound=1) + taper_start = validate_surface_threshold(interface_taper_start, + "interface_taper_start"; + upper_bound=1, + strict_upper_bound=true) + interpolation_threshold = validate_surface_threshold(interpolation_surface_threshold, + "interpolation_surface_threshold"; + upper_bound=1) + parameters = promote(boundary_threshold, normal_threshold, density_threshold, + taper_start, interpolation_threshold) + return eltype(parameters) <: Integer ? float.(parameters) : parameters +end + +function ColorfieldSurfaceDetection(; kwargs...) + return ColorfieldSurfaceDetection(colorfield_surface_parameters(; kwargs...)...) +end + +function ColorfieldSurfaceNormal(; kwargs...) + return ColorfieldSurfaceNormal(colorfield_surface_parameters(; kwargs...)...) +end + +function ColorfieldSurfaceNormal(boundary_contact_threshold, interface_threshold, + ideal_density_threshold) + return ColorfieldSurfaceNormal(; boundary_contact_threshold, interface_threshold, ideal_density_threshold) end -function create_cache_surface_normal(surface_normal_method, ELTYPE, NDIMS, nparticles) +@inline computes_surface_normal(surface_method) = false +@inline computes_surface_normal(::AbstractSurfaceNormalMethod) = true + +@inline is_colorfield_surface_method(surface_method) = false +@inline is_colorfield_surface_method(::ColorfieldSurfaceMethod) = true + +@inline contributes_boundary_colorfield(system) = false +@inline contributes_boundary_colorfield(::AbstractBoundarySystem) = true + +@inline function default_surface_method(surface_tension, surface_method) + if isnothing(surface_method) && requires_surface_normal(surface_tension) + return ColorfieldSurfaceNormal() + end + + return surface_method +end + +function select_surface_method(surface_tension, surface_method, surface_normal_method) + if !isnothing(surface_method) && !isnothing(surface_normal_method) + throw(ArgumentError("`surface_method` and deprecated `surface_normal_method` cannot both be set")) + end + + if !isnothing(surface_normal_method) + Base.depwarn("`surface_normal_method` is deprecated; use `surface_method` instead", + :surface_normal_method) + surface_method = surface_normal_method + end + + surface_method = default_surface_method(surface_tension, surface_method) + if !(surface_method isa Union{Nothing, AbstractSurfaceMethod}) + throw(ArgumentError("`surface_method` must be an `AbstractSurfaceMethod` or `nothing`")) + end + if requires_surface_normal(surface_tension) && !computes_surface_normal(surface_method) + throw(ArgumentError("$(typeof(surface_tension)) requires a surface method that computes surface normals")) + end + + return surface_method +end + +@inline function cubic_smoothstep(value) + value <= zero(value) && return zero(value) + value >= one(value) && return one(value) + return value^2 * (3 - 2value) +end + +@inline function gradient_surface_activity(normal_norm, support_radius, + surface_method::ColorfieldSurfaceMethod) + threshold = surface_method.interface_threshold + dimensionless_norm = support_radius * normal_norm + if iszero(threshold) + return iszero(dimensionless_norm) ? zero(dimensionless_norm) : + one(dimensionless_norm) + end + + lower_bound = surface_method.interface_taper_start * threshold + transition_coordinate = (dimensionless_norm - lower_bound) / + (threshold - lower_bound) + return cubic_smoothstep(transition_coordinate) +end + +function create_cache_surface(surface_method, ELTYPE, NDIMS, nparticles) return (;) end -function create_cache_surface_normal(::ColorfieldSurfaceNormal, ELTYPE, NDIMS, nparticles) +function create_cache_surface(::ColorfieldSurfaceDetection, ELTYPE, NDIMS, nparticles) + surface_gradient = Array{ELTYPE, 2}(undef, NDIMS, nparticles) + surface_activity = Array{ELTYPE, 1}(undef, nparticles) + neighbor_count = Array{ELTYPE, 1}(undef, nparticles) + colorfield = Array{ELTYPE, 1}(undef, nparticles) + return (; surface_gradient, surface_activity, neighbor_count, colorfield) +end + +function create_cache_surface(::ColorfieldSurfaceNormal, ELTYPE, NDIMS, nparticles) surface_normal = Array{ELTYPE, 2}(undef, NDIMS, nparticles) + surface_activity = Array{ELTYPE, 1}(undef, nparticles) neighbor_count = Array{ELTYPE, 1}(undef, nparticles) colorfield = Array{ELTYPE, 1}(undef, nparticles) correction_factor = Array{ELTYPE, 1}(undef, nparticles) - return (; surface_normal, neighbor_count, colorfield, correction_factor) + return (; surface_normal, surface_activity, neighbor_count, colorfield, + correction_factor) +end + +@inline function surface_gradient(cache, ::ColorfieldSurfaceDetection) + return cache.surface_gradient +end + +@inline function surface_gradient(cache, ::ColorfieldSurfaceNormal) + return cache.surface_normal end @inline function surface_normal(particle_system::AbstractFluidSystem, particle) - (; cache) = particle_system - return extract_svector(cache.surface_normal, particle_system, particle) + return extract_svector(particle_system.cache.surface_normal, particle_system, particle) +end + +@inline function surface_activity(particle_system::AbstractFluidSystem, particle) + return @inbounds particle_system.cache.surface_activity[particle] end -function calc_normal!(system, neighbor_system, u_system, v, v_neighbor_system, - u_neighbor_system, semi, surface_normal_method, - neighbor_surface_normal_method) - # Normal not needed +function calc_surface!(system, neighbor_system, u_system, v, v_neighbor_system, + u_neighbor_system, semi, surface_method, neighbor_surface_method) return system end # Section 2.2 in Akinci et al. 2013 "Versatile Surface Tension and Adhesion for SPH Fluids" # 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) - (; cache) = system +function calc_surface!(system::AbstractFluidSystem, + neighbor_system::AbstractFluidSystem, + u_system, v, v_neighbor_system, u_neighbor_system, semi, + surface_method::ColorfieldSurfaceMethod, + neighbor_surface_method) + contributes_to_colorfield(neighbor_system) || return system + (; cache) = system + gradient = surface_gradient(cache, surface_method) + color_b = neighbor_system.cache.color system_coords = current_coordinates(u_system, system) neighbor_system_coords = current_coordinates(u_neighbor_system, neighbor_system) @@ -60,61 +237,53 @@ function calc_normal!(system::AbstractFluidSystem, neighbor_system::AbstractFlui system_coords, neighbor_system_coords, semi; points=each_integrated_particle(system)) do particle, neighbor, pos_diff, distance - m_b = hydrodynamic_mass(neighbor_system, neighbor) - density_neighbor = current_density(v_neighbor_system, - neighbor_system, neighbor) - grad_kernel = smoothing_kernel_grad(system, pos_diff, distance, particle) + volume_b = hydrodynamic_mass(neighbor_system, neighbor) / + current_density(v_neighbor_system, neighbor_system, neighbor) + grad_kernel = kernel_grad(system_smoothing_kernel(system), pos_diff, distance, + smoothing_length(system, particle)) for i in 1:ndims(system) - cache.surface_normal[i, particle] += m_b / density_neighbor * grad_kernel[i] + gradient[i, particle] += volume_b * color_b * grad_kernel[i] end - cache.neighbor_count[particle] += 1 end return system end -# Section 2.2 in Akinci et al. 2013 "Versatile Surface Tension and Adhesion for SPH Fluids" -# Note: This is the simplest form of normal approximation commonly used in SPH and comes -# with serious deficits in accuracy especially at corners, small neighborhoods and boundaries -function calc_boundary_normal!(system::AbstractFluidSystem, neighbor_system, u_system, v, - u_neighbor_system, semi, surface_normal_method) +function calc_boundary_surface!(system::AbstractFluidSystem, neighbor_system, u_system, v, + u_neighbor_system, semi, + surface_method::ColorfieldSurfaceMethod) (; cache) = system + gradient = surface_gradient(cache, surface_method) (; colorfield, initial_colorfield) = neighbor_system.boundary_model.cache - (; boundary_contact_threshold) = surface_normal_method + (; boundary_contact_threshold) = surface_method + color_a = system.cache.color system_coords = current_coordinates(u_system, system) neighbor_system_coords = current_coordinates(u_neighbor_system, neighbor_system) - # First we need to calculate the smoothed colorfield values of the boundary - # 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 - - # Accumulate fluid neighbors + colorfield .= abs.(initial_colorfield) 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(color_a) * 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, semi) do particle, neighbor, pos_diff, distance - # We assume that we are in contact with the boundary if the color of the boundary particle - # is larger than the threshold if colorfield[neighbor] / maximum_colorfield > boundary_contact_threshold - m_b = hydrodynamic_mass(system, particle) - density_neighbor = current_density(v, system, particle) - grad_kernel = smoothing_kernel_grad(system, pos_diff, distance, particle) + volume_a = hydrodynamic_mass(system, particle) / + current_density(v, system, particle) + grad_kernel = kernel_grad(system_smoothing_kernel(system), pos_diff, distance, + smoothing_length(system, particle)) for i in 1:ndims(system) - cache.surface_normal[i, particle] += m_b / density_neighbor * grad_kernel[i] + gradient[i, particle] += volume_a * color_a * grad_kernel[i] end cache.neighbor_count[particle] += 1 end @@ -123,114 +292,128 @@ function calc_boundary_normal!(system::AbstractFluidSystem, neighbor_system, u_s return system end -function calc_normal!(system::AbstractFluidSystem, neighbor_system::AbstractBoundarySystem, - u_system, v, v_neighbor_system, u_neighbor_system, semi, - surface_normal_method, neighbor_surface_normal_method) - return calc_boundary_normal!(system, neighbor_system, u_system, v, u_neighbor_system, - semi, surface_normal_method) +function calc_surface!(system::AbstractFluidSystem, + neighbor_system::AbstractBoundarySystem, + u_system, v, v_neighbor_system, u_neighbor_system, semi, + surface_method::ColorfieldSurfaceMethod, + neighbor_surface_method) + return calc_boundary_surface!(system, neighbor_system, u_system, v, u_neighbor_system, + semi, surface_method) end -function remove_invalid_normals!(system::AbstractFluidSystem, surface_tension, - surface_normal_method) - (; cache) = system +function invalid_surface_particle(system, surface_method::ColorfieldSurfaceMethod, + particle, support_radius) + neighbor_count = system.cache.neighbor_count[particle] + minimum_neighbor_count = 2^ndims(system) + 1 + neighbor_count < minimum_neighbor_count && return true + + threshold = surface_method.ideal_density_threshold + return threshold > 0 && + threshold * ideal_neighbor_count(Val(ndims(system)), + system.cache.reference_particle_spacing, + support_radius) < neighbor_count +end + +function finalize_surface!(system::AbstractFluidSystem, surface_tension, + surface_method::ColorfieldSurfaceMethod) + gradient = surface_gradient(system.cache, surface_method) + support_radius = compact_support(system_smoothing_kernel(system), + initial_smoothing_length(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 + particle_gradient = extract_svector(gradient, system, particle) + gradient_norm = norm(particle_gradient) + activity = gradient_surface_activity(gradient_norm, support_radius, surface_method) + + if invalid_surface_particle(system, surface_method, particle, support_radius) + system.cache.surface_activity[particle] = zero(activity) + gradient[1:ndims(system), particle] .= 0 + else + system.cache.surface_activity[particle] = activity 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_normal_method::ColorfieldSurfaceNormal) - (; cache, smoothing_kernel) = system - (; ideal_density_threshold, interface_threshold) = surface_normal_method - (; neighbor_count) = cache - - smoothing_length_ = initial_smoothing_length(system) - - # We remove invalid normals i.e. they have a small norm (eq. 20) - normal_condition2 = (interface_threshold / - compact_support(smoothing_kernel, smoothing_length_))^2 +function finalize_surface!(system::AbstractFluidSystem, surface_tension, + surface_method::ColorfieldSurfaceNormal) + gradient = surface_gradient(system.cache, surface_method) + support_radius = compact_support(system_smoothing_kernel(system), + initial_smoothing_length(system)) + normal_condition2 = (surface_method.interface_threshold / support_radius)^2 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] - cache.surface_normal[1:ndims(system), particle] .= 0 - continue - end - - particle_surface_normal = surface_normal(system, particle) - norm2 = dot(particle_surface_normal, particle_surface_normal) - - # See eq. 21 - if norm2 > normal_condition2 - cache.surface_normal[1:ndims(system), - particle] = particle_surface_normal / sqrt(norm2) + particle_gradient = extract_svector(gradient, system, particle) + norm2 = dot(particle_gradient, particle_gradient) + gradient_norm = sqrt(norm2) + activity = gradient_surface_activity(gradient_norm, support_radius, surface_method) + + if invalid_surface_particle(system, surface_method, particle, support_radius) + system.cache.surface_activity[particle] = zero(activity) + gradient[1:ndims(system), particle] .= 0 + elseif norm2 > normal_condition2 + system.cache.surface_activity[particle] = activity + if normalize_surface_normals(surface_tension) + gradient[1:ndims(system), particle] = particle_gradient / gradient_norm + end else - cache.surface_normal[1:ndims(system), particle] .= 0 + system.cache.surface_activity[particle] = activity + gradient[1:ndims(system), particle] .= 0 end end return system end -function compute_surface_normal!(system, surface_normal_method, v, u, v_ode, u_ode, semi, t) +@inline normalize_surface_normals(surface_tension) = false +@inline normalize_surface_normals(::SurfaceTensionMorris) = true +@inline normalize_surface_normals(::SurfaceTensionMomentumMorris) = true + +function compute_surface!(system, surface_method, v, u, v_ode, u_ode, semi, t) return system end -function compute_surface_normal!(system::AbstractFluidSystem, - surface_normal_method_::ColorfieldSurfaceNormal, - v, u, v_ode, u_ode, semi, t) +function compute_surface!(system::AbstractFluidSystem, + surface_method_::ColorfieldSurfaceMethod, + v, u, v_ode, u_ode, semi, t) (; cache, surface_tension) = system - # Reset surface normal - set_zero!(cache.surface_normal) + set_zero!(surface_gradient(cache, surface_method_)) + set_zero!(cache.surface_activity) 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 + @trixi_timeit timer() "compute surface" begin foreach_system_wrapped(semi, v_ode, u_ode) do neighbor_system, v_neighbor_system, u_neighbor_system - if !has_system_interaction(system, neighbor_system, semi) - # No interaction between these systems. - return - end + has_system_interaction(system, neighbor_system, semi) || return - calc_normal!(system, neighbor_system, u, v, v_neighbor_system, - u_neighbor_system, semi, surface_normal_method_, - surface_normal_method(neighbor_system)) + calc_surface!(system, neighbor_system, u, v, v_neighbor_system, + u_neighbor_system, semi, surface_method_, + surface_method(neighbor_system)) end end - remove_invalid_normals!(system, surface_tension, surface_normal_method_) + finalize_surface!(system, surface_tension, surface_method_) return system end +function remove_invalid_normals!(system::AbstractFluidSystem, surface_tension, + surface_method::ColorfieldSurfaceNormal) + return finalize_surface!(system, surface_tension, surface_method) +end + function calc_curvature!(system, neighbor_system, u_system, v, - v_neighbor_system, u_neighbor_system, semi, surface_normal_method, - neighbor_surface_normal_method) + v_neighbor_system, u_neighbor_system, semi, surface_method, + neighbor_surface_method) end # Section 5 in Morris 2000 "Simulating surface tension with smoothed particle hydrodynamics" function calc_curvature!(system::AbstractFluidSystem, neighbor_system::AbstractFluidSystem, u_system, v, v_neighbor_system, u_neighbor_system, semi, - surface_normal_method::ColorfieldSurfaceNormal, - neighbor_surface_normal_method::ColorfieldSurfaceNormal) + surface_method_::ColorfieldSurfaceNormal, + neighbor_surface_method::ColorfieldSurfaceNormal) (; cache) = system (; curvature, correction_factor) = cache @@ -248,7 +431,6 @@ function calc_curvature!(system::AbstractFluidSystem, neighbor_system::AbstractF n_b = surface_normal(neighbor_system, neighbor) v_b = m_b / rho_b - # Eq. 22: we can test against `eps()` here since the surface normals that are invalid have been removed if dot(n_a, n_a) > eps() && dot(n_b, n_b) > eps() w = smoothing_kernel(system, distance, particle) grad_kernel = smoothing_kernel_grad(system, pos_diff, distance, particle) @@ -256,12 +438,10 @@ function calc_curvature!(system::AbstractFluidSystem, neighbor_system::AbstractF for i in 1:ndims(system) curvature[particle] += v_b * (n_b[i] - n_a[i]) * grad_kernel[i] end - # Eq. 24 correction_factor[particle] += v_b * w end end - # Eq. 23 for particle in each_integrated_particle(system) curvature[particle] /= (correction_factor[particle] + eps()) end @@ -276,23 +456,19 @@ end function compute_curvature!(system::AbstractFluidSystem, surface_tension::SurfaceTensionMorris, v, u, v_ode, u_ode, semi, t) - (; cache, surface_tension) = system + (; cache) = system - # Reset surface curvature set_zero!(cache.curvature) @trixi_timeit timer() "compute surface curvature" begin foreach_system_wrapped(semi, v_ode, u_ode) do neighbor_system, v_neighbor_system, u_neighbor_system - if !has_system_interaction(system, neighbor_system, semi) - # No interaction between these systems. - return - end + has_system_interaction(system, neighbor_system, semi) || return calc_curvature!(system, neighbor_system, u, v, v_neighbor_system, - u_neighbor_system, semi, surface_normal_method(system), - surface_normal_method(neighbor_system)) + u_neighbor_system, semi, surface_method(system), + surface_method(neighbor_system)) end end return system diff --git a/src/schemes/fluid/surface_tension.jl b/src/schemes/fluid/surface_tension.jl index 5656e95e12..3e0682328d 100644 --- a/src/schemes/fluid/surface_tension.jl +++ b/src/schemes/fluid/surface_tension.jl @@ -98,6 +98,11 @@ struct SurfaceTensionMomentumMorris{ELTYPE} <: AbstractSurfaceTension end end +@inline requires_surface_normal(surface_tension) = false +@inline requires_surface_normal(::SurfaceTensionAkinci) = true +@inline requires_surface_normal(::SurfaceTensionMorris) = true +@inline requires_surface_normal(::SurfaceTensionMomentumMorris) = true + 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 eea0607d7d..195325b61e 100644 --- a/src/schemes/fluid/weakly_compressible_sph/system.jl +++ b/src/schemes/fluid/weakly_compressible_sph/system.jl @@ -7,7 +7,7 @@ shifting_technique=nothing, buffer_size=nothing, correction=nothing, source_terms=nothing, - surface_tension=nothing, surface_normal_method=nothing, + surface_tension=nothing, surface_method=nothing, reference_particle_spacing=0.0, color_value=1)) System for particles of a fluid. @@ -53,14 +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) -- `reference_particle_spacing`: The reference particle spacing used for weighting values at the boundary, - which currently is only needed when using surface tension. -- `color_value`: Integer label used for calculation of surface normals. - Currently this is only used together with [`BoundaryModelDummyParticles`](@ref) and - [`ColorfieldSurfaceNormal`](@ref): fluid-boundary normal evaluation - reads the resulting boundary colorfield to detect wall contact. +- `surface_method`: Surface detection or normal method used by this system. + Methods that compute normals always also compute surface + activity. The default is `nothing`, or + `ColorfieldSurfaceNormal()` when required by surface tension. +- `reference_particle_spacing`: Reference spacing required by colorfield surface methods. +- `color_value`: Scalar contributed to colorfield surface detection. Different + values identify represented fluid-fluid interfaces. """ struct WeaklyCompressibleSPHSystem{NDIMS, ELTYPE <: Real, IC, MA, P, DC, SE, K, V, DD, COR, PF, SC, ST, B, SRFT, SRFN, PR, @@ -79,7 +78,7 @@ struct WeaklyCompressibleSPHSystem{NDIMS, ELTYPE <: Real, IC, MA, P, DC, SE, K, shifting_technique :: SC source_terms :: ST surface_tension :: SRFT - surface_normal_method :: SRFN + surface_method :: SRFN buffer :: B particle_refinement :: PR # TODO cache :: C @@ -97,7 +96,8 @@ function WeaklyCompressibleSPHSystem(initial_condition; smoothing_kernel, shifting_technique=nothing, buffer_size=nothing, correction=nothing, source_terms=nothing, - surface_tension=nothing, surface_normal_method=nothing, + surface_tension=nothing, surface_method=nothing, + surface_normal_method=nothing, reference_particle_spacing=0, color_value=1) buffer = isnothing(buffer_size) ? nothing : SystemBuffer(nparticles(initial_condition), buffer_size) @@ -130,12 +130,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_method = select_surface_method(surface_tension, surface_method, + 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")) + if is_colorfield_surface_method(surface_method) && reference_particle_spacing < eps() + throw(ArgumentError("`reference_particle_spacing` must be set to a positive value when using a colorfield surface method")) end pressure_acceleration = choose_pressure_acceleration_formulation(pressure_acceleration, @@ -146,8 +145,7 @@ function WeaklyCompressibleSPHSystem(initial_condition; smoothing_kernel, cache = (; create_cache_density(initial_condition, density_calculator)..., create_cache_correction(correction, initial_condition.density, NDIMS, n_particles)..., - create_cache_surface_normal(surface_normal_method, ELTYPE, NDIMS, - n_particles)..., + create_cache_surface(surface_method, ELTYPE, NDIMS, n_particles)..., create_cache_surface_tension(surface_tension, ELTYPE, NDIMS, n_particles)..., create_cache_refinement(initial_condition, particle_refinement, @@ -169,7 +167,7 @@ function WeaklyCompressibleSPHSystem(initial_condition; smoothing_kernel, smoothing_kernel, acceleration_, viscosity, density_diffusion, correction, pressure_acceleration, shifting_technique, source_terms, surface_tension, - surface_normal_method, buffer, particle_refinement, + surface_method, buffer, particle_refinement, cache) end @@ -185,8 +183,8 @@ function Base.show(io::IO, system::WeaklyCompressibleSPHSystem) print(io, ", ", system.density_diffusion) print(io, ", ", system.shifting_technique) print(io, ", ", system.surface_tension) - print(io, ", ", system.surface_normal_method) - if system.surface_normal_method isa ColorfieldSurfaceNormal + print(io, ", ", system.surface_method) + if is_colorfield_surface_method(system.surface_method) print(io, ", ", system.cache.color) end print(io, ", ", system.acceleration) @@ -217,8 +215,8 @@ function Base.show(io::IO, ::MIME"text/plain", system::WeaklyCompressibleSPHSyst summary_line(io, "density diffusion", system.density_diffusion) summary_line(io, "shifting technique", system.shifting_technique) summary_line(io, "surface tension", system.surface_tension) - summary_line(io, "surface normal method", system.surface_normal_method) - if system.surface_normal_method isa ColorfieldSurfaceNormal + summary_line(io, "surface method", system.surface_method) + if is_colorfield_surface_method(system.surface_method) summary_line(io, "color", system.cache.color) end summary_line(io, "acceleration", system.acceleration) @@ -321,7 +319,7 @@ end end function update_pressure!(system::WeaklyCompressibleSPHSystem, v, u, v_ode, u_ode, semi, t) - (; density_calculator, correction, surface_normal_method, surface_tension) = system + (; density_calculator, correction, surface_method, surface_tension) = system compute_pressure!(system, v, semi) @@ -332,8 +330,7 @@ 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 - compute_surface_normal!(system, surface_normal_method, v, u, v_ode, u_ode, semi, t) + compute_surface!(system, surface_method, v, u, v_ode, u_ode, semi, t) compute_surface_delta_function!(system, surface_tension, semi) return system end diff --git a/src/schemes/structure/rigid_body/system.jl b/src/schemes/structure/rigid_body/system.jl index dccacc2668..7ffb62ae59 100644 --- a/src/schemes/structure/rigid_body/system.jl +++ b/src/schemes/structure/rigid_body/system.jl @@ -32,9 +32,9 @@ torque and applied consistently to all rigid particles. 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 - colorfield initialization so fluids using - [`ColorfieldSurfaceNormal`](@ref) can detect contact with rigid + Currently this is used with `BoundaryModelDummyParticles` during + colorfield initialization so fluids using a colorfield surface method + 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"`. """ @@ -275,14 +275,18 @@ function initialize!(system::RigidBodySystem, semi) return system end -function calc_normal!(system::AbstractFluidSystem, - neighbor_system::RigidBodySystem{<:BoundaryModelDummyParticles}, - u_system, v, v_neighbor_system, u_neighbor_system, semi, - surface_normal_method, neighbor_surface_normal_method) +function calc_surface!(system::AbstractFluidSystem, + neighbor_system::RigidBodySystem{<:BoundaryModelDummyParticles}, + u_system, v, v_neighbor_system, u_neighbor_system, semi, + surface_method_, neighbor_surface_method) haskey(neighbor_system.boundary_model.cache, :initial_colorfield) || return system - return calc_boundary_normal!(system, neighbor_system, u_system, v, u_neighbor_system, - semi, surface_normal_method) + return calc_boundary_surface!(system, neighbor_system, u_system, v, u_neighbor_system, + semi, surface_method_) +end + +@inline function contributes_boundary_colorfield(::RigidBodySystem{<:BoundaryModelDummyParticles}) + return true end @inline function adhesion_force!(dv_particle, @@ -631,10 +635,10 @@ function check_configuration(system::RigidBodySystem, systems, nhs) end if neighbor isa AbstractFluidSystem && - neighbor.surface_normal_method isa ColorfieldSurfaceNormal + is_colorfield_surface_method(surface_method(neighbor)) if !(boundary_model isa BoundaryModelDummyParticles) throw(ArgumentError("`RigidBodySystem` is only compatible with " * - "`ColorfieldSurfaceNormal` when using " * + "colorfield surface methods when using " * "`BoundaryModelDummyParticles`.")) end @@ -642,7 +646,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.")) + "a colorfield surface method or a surface tension model.")) end end end diff --git a/test/general/interpolation.jl b/test/general/interpolation.jl index 5469876724..bd8420ecd4 100644 --- a/test/general/interpolation.jl +++ b/test/general/interpolation.jl @@ -131,6 +131,111 @@ semi_boundary = Semidiscretization(fluid_system, boundary_system) TrixiParticles.initialize_neighborhood_searches!(semi_boundary) + surface_detection_ic = RectangularShape(particle_spacing, (nx, ny), (0.0, 0.0), + density=1000.0) + surface_detection_system = WeaklyCompressibleSPHSystem(surface_detection_ic; + smoothing_kernel, + smoothing_length=1.5 * + particle_spacing, + density_calculator=ContinuityDensity(), + state_equation, viscosity, + acceleration=(0.0, -9.81), + surface_method=ColorfieldSurfaceDetection(ideal_density_threshold=0.9), + reference_particle_spacing=particle_spacing) + surface_detection_system.pressure .= surface_detection_ic.pressure + semi_surface_detection = Semidiscretization(surface_detection_system) + TrixiParticles.initialize_neighborhood_searches!(semi_surface_detection) + + @testset verbose=true "Interpolated Free Surface Detection" begin + min_x, max_x = extrema(view(surface_detection_ic.coordinates, 1, :)) + max_y = maximum(view(surface_detection_ic.coordinates, 2, :)) + point_coords = [(min_x + max_x) / 2 (min_x + max_x) / 2; + max_y max_y + particle_spacing] + v_surface_detection = vcat(surface_detection_ic.velocity, + surface_detection_ic.density') + u_surface_detection = surface_detection_ic.coordinates + + uncut = interpolate_points(point_coords, semi_surface_detection, + surface_detection_system, v_surface_detection, + u_surface_detection; + cut_off_bnd=false) + cut = interpolate_points(point_coords, semi_surface_detection, + surface_detection_system, v_surface_detection, + u_surface_detection; + cut_off_bnd=true) + + @test all(isfinite, uncut.density) + @test all(isfinite, uncut.surface_activity) + @test uncut.surface_activity[2] > 0.9 + @test isfinite(cut.density[1]) + @test cut.neighbor_count[1] > 0 + @test isnan(cut.density[2]) + @test isnan(cut.surface_activity[2]) + @test cut.neighbor_count[2] == 0 + + mktempdir() do output_directory + interpolate_plane_2d_vtk([min_x, max_y - particle_spacing], + [max_x, max_y + particle_spacing], + particle_spacing, semi_surface_detection, + surface_detection_system, v_surface_detection, + u_surface_detection; output_directory, + filename="surface_detection") + vtk_file = TrixiParticles.ReadVTK.VTKFile(joinpath(output_directory, + "surface_detection.vti")) + point_data = TrixiParticles.ReadVTK.get_point_data(vtk_file) + @test "surface_activity" in keys(point_data) + end + end + + @testset verbose=true "Interpolated Multiphase Surface Detection" begin + interface_spacing = 0.1 + y_coordinates = collect(-0.5:interface_spacing:0.5) + coordinates_a = hcat(([x, y] for x in -0.5:interface_spacing:-0.1 + for y in y_coordinates)...) + coordinates_b = hcat(([x, y] for x in 0.0:interface_spacing:0.5 + for y in y_coordinates)...) + initial_condition_a = InitialCondition(; coordinates=coordinates_a, + density=fill(1000.0, + size(coordinates_a, 2)), + particle_spacing=interface_spacing) + initial_condition_b = InitialCondition(; coordinates=coordinates_b, + density=fill(1000.0, + size(coordinates_b, 2)), + particle_spacing=interface_spacing) + interface_kernel = WendlandC2Kernel{2}() + interface_state_equation = StateEquationCole(sound_speed=10.0, + reference_density=1000.0, + exponent=1) + system_a = WeaklyCompressibleSPHSystem(initial_condition_a; + smoothing_kernel=interface_kernel, + smoothing_length=0.15, + density_calculator=SummationDensity(), + state_equation=interface_state_equation, + surface_method=ColorfieldSurfaceDetection(interface_threshold=1.0e-6), + reference_particle_spacing=interface_spacing, + color_value=1) + system_b = WeaklyCompressibleSPHSystem(initial_condition_b; + smoothing_kernel=interface_kernel, + smoothing_length=0.15, + density_calculator=SummationDensity(), + state_equation=interface_state_equation, + color_value=2) + interface_semi = Semidiscretization(system_a, system_b) + interface_ode = semidiscretize(interface_semi, (0.0, 0.01)) + TrixiParticles.update_systems_and_nhs(interface_ode.u0.x..., interface_semi, + 0.0) + interface_points = [-0.1 0.1; 0.0 0.0] + uncut = interpolate_points(interface_points, interface_semi, system_a, + interface_ode.u0.x...; cut_off_bnd=false) + cut = interpolate_points(interface_points, interface_semi, system_a, + interface_ode.u0.x...; cut_off_bnd=true) + + @test all(>(0.9), uncut.surface_activity) + @test isfinite(cut.density[1]) + @test isnan(cut.density[2]) + @test isnan(cut.surface_activity[2]) + end + # Some simple results expected_zero(y) = (density=[NaN], neighbor_count=[0], point_coords=[0.0; y;;], velocity=[NaN; NaN;;], pressure=[NaN]) diff --git a/test/general/semidiscretization.jl b/test/general/semidiscretization.jl index 5b92d133ba..b7a0e603a4 100644 --- a/test/general/semidiscretization.jl +++ b/test/general/semidiscretization.jl @@ -132,7 +132,7 @@ # Mock fluid system struct FluidSystemMock <: TrixiParticles.AbstractFluidSystem{2} surface_tension::Nothing - surface_normal_method::Nothing + surface_method::Nothing FluidSystemMock() = new(nothing, nothing) end diff --git a/test/schemes/fluid/surface_normal_sph.jl b/test/schemes/fluid/surface_normal_sph.jl index 5eb8a81704..cb0e770051 100644 --- a/test/schemes/fluid/surface_normal_sph.jl +++ b/test/schemes/fluid/surface_normal_sph.jl @@ -61,7 +61,7 @@ end function create_fluid_system(coordinates, velocity, mass, density, particle_spacing, surface_tension; - surface_normal_method=ColorfieldSurfaceNormal(), + surface_method=ColorfieldSurfaceNormal(), color_value=1, NDIMS=2, smoothing_length=1.0, wall=false, walldistance=0.0, boundary_system_type=:wall, smoothing_kernel=SchoenbergCubicSplineKernel{NDIMS}()) @@ -76,9 +76,9 @@ function create_fluid_system(coordinates, velocity, mass, density, particle_spac system = WeaklyCompressibleSPHSystem(fluid; smoothing_kernel, smoothing_length, density_calculator=SummationDensity(), state_equation, - surface_normal_method, + surface_method, reference_particle_spacing=particle_spacing, - surface_tension) + surface_tension, color_value) if wall boundary_system = if boundary_system_type == :wall @@ -109,12 +109,11 @@ function compute_and_test_surface_values(system, semi, ode; NDIMS=2) v = TrixiParticles.wrap_v(v0_ode, system, semi) u = TrixiParticles.wrap_u(u0_ode, system, semi) - # Compute the surface normals - TrixiParticles.compute_surface_normal!(system, system.surface_normal_method, v, u, - v0_ode, u0_ode, semi, 0.0) + TrixiParticles.compute_surface!(system, system.surface_method, v, u, + v0_ode, u0_ode, semi, 0.0) TrixiParticles.remove_invalid_normals!(system, system.surface_tension, - system.surface_normal_method) + system.surface_method) # After computation, check that surface normals have been computed and are not NaN or Inf @test all(isfinite, system.cache.surface_normal) @@ -143,6 +142,225 @@ function compute_curvature!(system, semi, ode) v, u, v0_ode, u0_ode, semi, 0.0) end +@testset verbose=true "Colorfield Surface Detection" begin + normal_method = ColorfieldSurfaceNormal(ideal_density_threshold=0.9) + detection_method = ColorfieldSurfaceDetection(ideal_density_threshold=0.9) + @test ColorfieldSurfaceNormal() == ColorfieldSurfaceNormal(0.1, 0.01, 0.0) + @test normal_method.interpolation_surface_threshold == 0.45 + @test detection_method.interpolation_surface_threshold == 0.45 + @test TrixiParticles.computes_surface_normal(normal_method) + @test !TrixiParticles.computes_surface_normal(detection_method) + + @test_throws ArgumentError ColorfieldSurfaceNormal(boundary_contact_threshold=-0.1) + @test_throws ArgumentError ColorfieldSurfaceNormal(interface_threshold=Inf) + @test_throws ArgumentError ColorfieldSurfaceNormal(interface_taper_start=1.0) + @test_throws ArgumentError ColorfieldSurfaceNormal(interpolation_surface_threshold=1.1) + @test_throws ArgumentError ColorfieldSurfaceDetection(interface_threshold=-0.1) + @test_throws ArgumentError ColorfieldSurfaceDetection(interface_threshold="invalid") + + particle_spacing = 0.1 + coordinates = RectangularShape(particle_spacing, (21, 11), (0.0, 0.0), + density=1000.0) + smoothing_kernel = SchoenbergCubicSplineKernel{2}() + system, _, semi, + ode = create_fluid_system(coordinates.coordinates, coordinates.velocity, + coordinates.mass, coordinates.density, particle_spacing, + nothing; smoothing_length=1.5 * particle_spacing, + smoothing_kernel, surface_method=normal_method) + + activity = system.cache.surface_activity + x = coordinates.coordinates + min_x, max_x = extrema(view(x, 1, :)) + min_y, max_y = extrema(view(x, 2, :)) + interior = [particle + for particle in eachindex(activity) + if min_x + 3particle_spacing < x[1, particle] < + max_x - 3particle_spacing && + min_y + 3particle_spacing < x[2, particle] < + max_y - 3particle_spacing] + top_surface = [particle + for particle in eachindex(activity) + if x[2, particle] == max_y && + min_x + 3particle_spacing < x[1, particle] < + max_x - 3particle_spacing] + + @test all(isfinite, activity) + @test all(iszero, activity[interior]) + @test all(==(1), activity[top_surface]) + + detection_system = WeaklyCompressibleSPHSystem(coordinates; + smoothing_kernel, + smoothing_length=1.5 * particle_spacing, + density_calculator=SummationDensity(), + state_equation=system.state_equation, + surface_method=detection_method, + reference_particle_spacing=particle_spacing) + detection_semi = Semidiscretization(detection_system) + detection_ode = semidiscretize(detection_semi, (0.0, 0.01)) + TrixiParticles.update_systems_and_nhs(detection_ode.u0.x..., detection_semi, 0.0) + @test !haskey(detection_system.cache, :surface_normal) + @test detection_system.cache.surface_activity == activity + @test isapprox(detection_system.cache.surface_gradient, system.cache.surface_normal; + rtol=10eps(), atol=10eps()) + + edac_system = EntropicallyDampedSPHSystem(coordinates; + smoothing_kernel, + smoothing_length=1.5 * particle_spacing, + sound_speed=10.0, + density_calculator=SummationDensity(), + surface_method=detection_method, + reference_particle_spacing=particle_spacing) + edac_semi = Semidiscretization(edac_system) + edac_ode = semidiscretize(edac_semi, (0.0, 0.01)) + TrixiParticles.update_systems_and_nhs(edac_ode.u0.x..., edac_semi, 0.0) + @test edac_system.cache.surface_activity == activity + + iisph_system = ImplicitIncompressibleSPHSystem(coordinates; + smoothing_kernel, + smoothing_length=1.5 * particle_spacing, + reference_density=1000.0, + time_step=0.001, + surface_method=detection_method, + reference_particle_spacing=particle_spacing) + iisph_semi = Semidiscretization(iisph_system) + iisph_ode = semidiscretize(iisph_semi, (0.0, 0.01)) + TrixiParticles.update_systems_and_nhs(iisph_ode.u0.x..., iisph_semi, 0.0) + @test iisph_system.cache.surface_activity == activity + + @test_throws ArgumentError WeaklyCompressibleSPHSystem(coordinates; + smoothing_kernel, + smoothing_length=1.5 * + particle_spacing, + density_calculator=SummationDensity(), + state_equation=system.state_equation, + surface_tension=SurfaceTensionMorris(), + surface_method=detection_method, + reference_particle_spacing=particle_spacing) + cohesion_system = WeaklyCompressibleSPHSystem(coordinates; + smoothing_kernel, + smoothing_length=1.5 * particle_spacing, + density_calculator=SummationDensity(), + state_equation=system.state_equation, + surface_tension=CohesionForceAkinci()) + @test isnothing(cohesion_system.surface_method) + @test !haskey(cohesion_system.cache, :surface_activity) + + corrected_system = WeaklyCompressibleSPHSystem(coordinates; + smoothing_kernel, + smoothing_length=1.5 * particle_spacing, + density_calculator=SummationDensity(), + state_equation=system.state_equation, + correction=GradientCorrection(), + surface_method=normal_method, + reference_particle_spacing=particle_spacing) + corrected_semi = Semidiscretization(corrected_system) + corrected_ode = semidiscretize(corrected_semi, (0.0, 0.01)) + TrixiParticles.update_systems_and_nhs(corrected_ode.u0.x..., corrected_semi, 0.0) + @test corrected_system.cache.surface_activity == activity + @test isapprox(corrected_system.cache.surface_normal, system.cache.surface_normal; + rtol=10eps(), atol=10eps()) + + v_ode, u_ode = ode.u0.x + @test TrixiParticles.surface_activity(system, nothing, nothing, v_ode, u_ode, + semi, 0.0) == activity + @test TrixiParticles.surface_normal(detection_system, nothing, nothing, + detection_ode.u0.x..., detection_semi, 0.0) === + nothing + + metadata = Dict{String, Any}() + TrixiParticles.add_system_data!(metadata, normal_method) + @test metadata["surface_method"]["computes_surface_normal"] + @test metadata["surface_method"]["interpolation_surface_threshold"] == 0.45 + + mktempdir() do output_directory + trixi2vtk(ode.u0, semi, 0.0; output_directory, + prefix="surface_detection", overwrite=true) + vtk_data = vtk2trixi(joinpath(output_directory, + "surface_detection_fluid_1_current.vtu")) + @test vtk_data.surface_activity == activity + @test hasproperty(vtk_data, :surf_normal) + + trixi2vtk(detection_ode.u0, detection_semi, 0.0; output_directory, + prefix="detection_only", overwrite=true) + detection_vtk_data = vtk2trixi(joinpath(output_directory, + "detection_only_fluid_1_current.vtu")) + @test detection_vtk_data.surface_activity == activity + @test !hasproperty(detection_vtk_data, :surf_normal) + end +end + +@testset verbose=true "Multicolor Surface Activity" 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) + + function interface_geometry(color_a, color_b, surface_method_) + 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_method=surface_method_, + reference_particle_spacing=particle_spacing, + color_value=color_a) + 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)) + TrixiParticles.update_systems_and_nhs(ode.u0.x..., 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 + gradient = if surface_method_ isa ColorfieldSurfaceNormal + TrixiParticles.surface_normal(system_a, interface_particle) + else + TrixiParticles.extract_svector(system_a.cache.surface_gradient, system_a, + interface_particle) + end + return gradient, TrixiParticles.surface_activity(system_a, interface_particle), + system_b + end + + detection_method = ColorfieldSurfaceDetection(interface_threshold=1.0e-6) + normal_method = ColorfieldSurfaceNormal(interface_threshold=1.0e-6) + increasing_gradient, increasing_activity, + non_surface_neighbor = interface_geometry(0, 2, detection_method) + unit_gradient, _, _ = interface_geometry(0, 1, detection_method) + decreasing_gradient, decreasing_activity, _ = interface_geometry(2, 0, + detection_method) + equal_gradient, equal_activity, _ = interface_geometry(1, 1, detection_method) + normal_gradient, normal_activity, _ = interface_geometry(0, 2, normal_method) + + @test isnothing(non_surface_neighbor.surface_method) + @test increasing_gradient[1] > 0 + @test decreasing_gradient[1] < 0 + @test isapprox(norm(increasing_gradient), 2norm(unit_gradient); rtol=1.0e-12) + @test norm(equal_gradient) < 100eps() + @test increasing_activity == 1 + @test decreasing_activity == 1 + @test equal_activity == 0 + @test normal_activity == increasing_activity + @test normal_gradient == increasing_gradient +end + @testset verbose=true "Rigid Dummy Boundary Matches Wall Boundary" begin NDIMS = 2 particle_spacing = 0.2 @@ -161,16 +379,16 @@ end wall_ode = create_fluid_system(coordinates, velocity, mass, density, particle_spacing, SurfaceTensionMorris(surface_tension_coefficient=0.072); NDIMS, smoothing_length, smoothing_kernel, - surface_normal_method=ColorfieldSurfaceNormal(interface_threshold=0.1, - ideal_density_threshold=0.9), + surface_method=ColorfieldSurfaceNormal(interface_threshold=0.1, + ideal_density_threshold=0.9), wall=true, walldistance=2.0, boundary_system_type=:wall) rigid_system, rigid_boundary, rigid_semi, rigid_ode = create_fluid_system(coordinates, velocity, mass, density, particle_spacing, SurfaceTensionMorris(surface_tension_coefficient=0.072); NDIMS, smoothing_length, smoothing_kernel, - surface_normal_method=ColorfieldSurfaceNormal(interface_threshold=0.1, - ideal_density_threshold=0.9), + surface_method=ColorfieldSurfaceNormal(interface_threshold=0.1, + ideal_density_threshold=0.9), wall=true, walldistance=2.0, boundary_system_type=:rigid) @@ -186,6 +404,9 @@ end @test isapprox(rigid_system.cache.neighbor_count, wall_system.cache.neighbor_count, rtol=sqrt(eps()), atol=sqrt(eps())) + @test isapprox(rigid_system.cache.surface_activity, + wall_system.cache.surface_activity, + rtol=sqrt(eps()), atol=sqrt(eps())) end @testset verbose=true "CSS/CSF: Sphere Surface Normals" begin @@ -219,8 +440,8 @@ end particle_spacing, SurfaceTensionMorris(surface_tension_coefficient=0.072); NDIMS, smoothing_length, smoothing_kernel, - surface_normal_method=ColorfieldSurfaceNormal(interface_threshold=0.1, - ideal_density_threshold=0.9), + surface_method=ColorfieldSurfaceNormal(interface_threshold=0.1, + ideal_density_threshold=0.9), wall=true, walldistance=2.0) compute_and_test_surface_values(system, semi, ode; NDIMS) @@ -319,8 +540,8 @@ end particle_spacing, SurfaceTensionAkinci(surface_tension_coefficient=0.072); NDIMS, smoothing_length, smoothing_kernel, - surface_normal_method=ColorfieldSurfaceNormal(interface_threshold=0.1, - ideal_density_threshold=0.9), + surface_method=ColorfieldSurfaceNormal(interface_threshold=0.1, + ideal_density_threshold=0.9), wall=true, walldistance=2.0) compute_and_test_surface_values(system, semi, ode; NDIMS) diff --git a/test/schemes/fluid/surface_tension.jl b/test/schemes/fluid/surface_tension.jl index 7fe8abbd97..404d934552 100644 --- a/test/schemes/fluid/surface_tension.jl +++ b/test/schemes/fluid/surface_tension.jl @@ -115,8 +115,8 @@ 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), + surface_method=ColorfieldSurfaceNormal(interface_threshold=0.1, + ideal_density_threshold=0.9), reference_particle_spacing=1.0,) # 4. Verify Cache Contains Necessary Fields diff --git a/test/systems/edac_system.jl b/test/systems/edac_system.jl index ce04cd6774..511f691351 100644 --- a/test/systems/edac_system.jl +++ b/test/systems/edac_system.jl @@ -143,7 +143,7 @@ │ average pressure reduction: ……… no │ │ acceleration: …………………………………………… [0.0, 0.0] │ │ surface tension: …………………………………… nothing │ - │ surface normal method: …………………… nothing │ + │ surface method: ……………………………………… nothing │ └──────────────────────────────────────────────────────────────────────────────────────────────────┘""" @test repr("text/plain", system) == show_box end diff --git a/test/systems/wcsph_system.jl b/test/systems/wcsph_system.jl index 37b94f3f64..1dece49a04 100644 --- a/test/systems/wcsph_system.jl +++ b/test/systems/wcsph_system.jl @@ -213,7 +213,7 @@ │ density diffusion: ……………………………… Val{:density_diffusion}() │ │ shifting technique: …………………………… nothing │ │ surface tension: …………………………………… nothing │ - │ surface normal method: …………………… nothing │ + │ surface method: ……………………………………… nothing │ │ acceleration: …………………………………………… [0.0, 0.0] │ │ source terms: …………………………………………… Nothing │ └──────────────────────────────────────────────────────────────────────────────────────────────────┘"""