From 393f9dc907ca822ef5b180889d5d92b076cb142d Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Sat, 16 Aug 2025 12:24:13 +0300 Subject: [PATCH 1/3] Make average pressure reduction for EDAC independent of TVF --- .../fluid/entropically_damped_sph/rhs.jl | 10 ++-- .../fluid/entropically_damped_sph/system.jl | 54 +++++++++++++------ src/schemes/fluid/transport_velocity.jl | 15 +----- .../fluid/weakly_compressible_sph/system.jl | 2 +- test/systems/edac_system.jl | 33 +++++++----- 5 files changed, 67 insertions(+), 47 deletions(-) diff --git a/src/schemes/fluid/entropically_damped_sph/rhs.jl b/src/schemes/fluid/entropically_damped_sph/rhs.jl index 7a0d1ded58..08139e2ca6 100644 --- a/src/schemes/fluid/entropically_damped_sph/rhs.jl +++ b/src/schemes/fluid/entropically_damped_sph/rhs.jl @@ -27,10 +27,12 @@ function interact!(dv, v_particle_system, u_particle_system, p_a = current_pressure(v_particle_system, particle_system, particle) p_b = current_pressure(v_neighbor_system, neighbor_system, neighbor) - # This technique is for a more robust `pressure_acceleration` but only with TVF. - # It results only in significant improvement for EDAC and not for WCSPH. - # See Ramachandran (2019) p. 582 - # Note that the return value is zero when not using EDAC with TVF. + # This technique by Basa et al. 2017 (10.1002/fld.1927) aims to reduce numerical + # errors due to large pressures by subtracting the average pressure of neighboring + # particles. + # It results in significant improvement for EDAC, especially with TVF, + # but not for WCSPH, according to Ramachandran & Puri (2019), Section 3.2. + # Note that the return value is zero when not using average pressure reduction. p_avg = average_pressure(particle_system, particle) m_a = hydrodynamic_mass(particle_system, particle) diff --git a/src/schemes/fluid/entropically_damped_sph/system.jl b/src/schemes/fluid/entropically_damped_sph/system.jl index f9875e684a..d197e5ef02 100644 --- a/src/schemes/fluid/entropically_damped_sph/system.jl +++ b/src/schemes/fluid/entropically_damped_sph/system.jl @@ -30,7 +30,10 @@ See [Entropically Damped Artificial Compressibility for SPH](@ref edac) for more When set to `nothing`, the pressure acceleration formulation for the corresponding [density calculator](@ref density_calculator) is chosen. - `density_calculator`: [Density calculator](@ref density_calculator) (default: [`SummationDensity`](@ref)) -- `transport_velocity`: [Transport Velocity Formulation (TVF)](@ref transport_velocity_formulation). Default is no TVF. +- `transport_velocity`: [Transport Velocity Formulation (TVF)](@ref transport_velocity_formulation). + Default is no TVF. +- `average_pressure_reduction`: Whether to subtract the average pressure of neighboring particles + from the local pressure (default: `true` when using TVF, `false` otherwise). - `buffer_size`: Number of buffer particles. This is needed when simulating with [`OpenBoundarySPHSystem`](@ref). - `correction`: Correction method used for this system. (default: no correction, see [Corrections](@ref corrections)) @@ -53,7 +56,7 @@ See [Entropically Damped Artificial Compressibility for SPH](@ref edac) for more """ struct EntropicallyDampedSPHSystem{NDIMS, ELTYPE <: Real, IC, M, DC, K, V, COR, PF, TV, - ST, SRFT, SRFN, B, PR, C} <: FluidSystem{NDIMS} + AVGP, ST, SRFT, SRFN, B, PR, C} <: FluidSystem{NDIMS} initial_condition :: IC mass :: M # Vector{ELTYPE}: [particle] density_calculator :: DC @@ -65,6 +68,7 @@ struct EntropicallyDampedSPHSystem{NDIMS, ELTYPE <: Real, IC, M, DC, K, V, COR, correction :: COR pressure_acceleration_formulation :: PF transport_velocity :: TV + average_pressure_reduction :: AVGP source_terms :: ST surface_tension :: SRFT surface_normal_method :: SRFN @@ -80,6 +84,7 @@ function EntropicallyDampedSPHSystem(initial_condition, smoothing_kernel, pressure_acceleration=inter_particle_averaged_pressure, density_calculator=SummationDensity(), transport_velocity=nothing, + average_pressure_reduction=!isnothing(transport_velocity), alpha=0.5, viscosity=nothing, acceleration=ntuple(_ -> 0.0, ndims(smoothing_kernel)), @@ -127,10 +132,14 @@ function EntropicallyDampedSPHSystem(initial_condition, smoothing_kernel, NDIMS, ELTYPE, correction) + avg_pressure_reduction = Val(average_pressure_reduction) + nu_edac = (alpha * smoothing_length * sound_speed) / 8 cache = (; create_cache_density(initial_condition, density_calculator)..., - create_cache_tvf_edac(initial_condition, transport_velocity)..., + create_cache_tvf(initial_condition, transport_velocity)..., + create_cache_avg_pressure_reduction(initial_condition, + avg_pressure_reduction)..., create_cache_surface_normal(surface_normal_method, ELTYPE, NDIMS, n_particles)..., create_cache_surface_tension(surface_tension, ELTYPE, NDIMS, @@ -152,19 +161,29 @@ function EntropicallyDampedSPHSystem(initial_condition, smoothing_kernel, EntropicallyDampedSPHSystem{NDIMS, ELTYPE, typeof(initial_condition), typeof(mass), typeof(density_calculator), typeof(smoothing_kernel), typeof(viscosity), typeof(correction), - typeof(pressure_acceleration), - typeof(transport_velocity), typeof(source_terms), + typeof(pressure_acceleration), typeof(transport_velocity), + typeof(avg_pressure_reduction), typeof(source_terms), typeof(surface_tension), typeof(surface_normal_method), typeof(buffer), Nothing, typeof(cache)}(initial_condition, mass, density_calculator, smoothing_kernel, sound_speed, viscosity, nu_edac, acceleration_, correction, pressure_acceleration, transport_velocity, + avg_pressure_reduction, source_terms, surface_tension, surface_normal_method, buffer, particle_refinement, cache) end +create_cache_avg_pressure_reduction(initial_condition, ::Val{false}) = (;) + +function create_cache_avg_pressure_reduction(initial_condition, ::Val{true}) + pressure_average = copy(initial_condition.pressure) + neighbor_counter = Vector{Int}(undef, nparticles(initial_condition)) + + return (; pressure_average, neighbor_counter) +end + function Base.show(io::IO, system::EntropicallyDampedSPHSystem) @nospecialize system # reduce precompilation time @@ -201,6 +220,8 @@ function Base.show(io::IO, ::MIME"text/plain", system::EntropicallyDampedSPHSyst summary_line(io, "smoothing kernel", system.smoothing_kernel |> typeof |> nameof) summary_line(io, "tansport velocity formulation", system.transport_velocity |> typeof |> nameof) + summary_line(io, "average pressure reduction", + 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) @@ -252,17 +273,15 @@ end @inline transport_velocity(system::EntropicallyDampedSPHSystem) = system.transport_velocity -@inline average_pressure(system, particle) = zero(eltype(system)) - @inline function average_pressure(system::EntropicallyDampedSPHSystem, particle) - average_pressure(system, transport_velocity(system), particle) + average_pressure(system, system.average_pressure_reduction, particle) end -@inline function average_pressure(system, ::TransportVelocityAdami, particle) +@inline function average_pressure(system, ::Val{true}, particle) return system.cache.pressure_average[particle] end -@inline average_pressure(system, ::Nothing, particle) = zero(eltype(system)) +@inline average_pressure(system, ::Val{false}, particle) = zero(eltype(system)) @inline function current_density(v, system::EntropicallyDampedSPHSystem) return current_density(v, system.density_calculator, system) @@ -301,18 +320,21 @@ function update_final!(system::EntropicallyDampedSPHSystem, v, u, v_ode, u_ode, # Surface normal of neighbor and boundary needs to have been calculated already compute_curvature!(system, surface_tension, v, u, v_ode, u_ode, semi, t) compute_stress_tensors!(system, surface_tension, v, u, v_ode, u_ode, semi, t) - update_average_pressure!(system, transport_velocity(system), v_ode, u_ode, semi) + update_average_pressure!(system, system.average_pressure_reduction, v_ode, u_ode, semi) update_tvf!(system, transport_velocity(system), v, u, v_ode, u_ode, semi, t) end -function update_average_pressure!(system, ::Nothing, v_ode, u_ode, semi) +# No average pressure reduction is used +function update_average_pressure!(system, ::Val{false}, v_ode, u_ode, semi) return system end -# This technique is for a more robust `pressure_acceleration` but only with TVF. -# It results only in significant improvement for EDAC and not for WCSPH. -# See Ramachandran (2019) p. 582. -function update_average_pressure!(system, ::TransportVelocityAdami, v_ode, u_ode, semi) +# This technique by Basa et al. 2017 (10.1002/fld.1927) aims to reduce numerical errors +# due to large pressures by subtracting the average pressure of neighboring particles. +# It results in significant improvement for EDAC, especially with TVF, +# but not for WCSPH, according to Ramachandran & Puri (2019), Section 3.2. +# See eq. 16 and 17 in Ramachandran & Puri (2019) for an explanation of the technique. +function update_average_pressure!(system, ::Val{true}, v_ode, u_ode, semi) (; cache) = system (; pressure_average, neighbor_counter) = cache diff --git a/src/schemes/fluid/transport_velocity.jl b/src/schemes/fluid/transport_velocity.jl index bbf40bde2c..fce7375dd5 100644 --- a/src/schemes/fluid/transport_velocity.jl +++ b/src/schemes/fluid/transport_velocity.jl @@ -16,26 +16,15 @@ end # No TVF for a system by default @inline transport_velocity(system) = nothing -create_cache_tvf_wcsph(initial_condition, ::Nothing) = (;) +create_cache_tvf(initial_condition, ::Nothing) = (;) -function create_cache_tvf_wcsph(initial_condition, ::TransportVelocityAdami) +function create_cache_tvf(initial_condition, ::TransportVelocityAdami) delta_v = zeros(eltype(initial_condition), ndims(initial_condition), nparticles(initial_condition)) return (; delta_v) end -create_cache_tvf_edac(initial_condition, ::Nothing) = (;) - -function create_cache_tvf_edac(initial_condition, ::TransportVelocityAdami) - delta_v = zeros(eltype(initial_condition), ndims(initial_condition), - nparticles(initial_condition)) - pressure_average = copy(initial_condition.pressure) - neighbor_counter = Vector{Int}(undef, nparticles(initial_condition)) - - return (; delta_v, pressure_average, neighbor_counter) -end - # `δv` is the correction to the particle velocity due to the TVF. # Particles are advected with the velocity `v + δv`. @propagate_inbounds function delta_v(system, particle) diff --git a/src/schemes/fluid/weakly_compressible_sph/system.jl b/src/schemes/fluid/weakly_compressible_sph/system.jl index bc9d541ea8..a5f56499b8 100644 --- a/src/schemes/fluid/weakly_compressible_sph/system.jl +++ b/src/schemes/fluid/weakly_compressible_sph/system.jl @@ -147,7 +147,7 @@ function WeaklyCompressibleSPHSystem(initial_condition, n_particles)..., create_cache_refinement(initial_condition, particle_refinement, smoothing_length)..., - create_cache_tvf_wcsph(initial_condition, transport_velocity)..., + create_cache_tvf(initial_condition, transport_velocity)..., color=Int(color_value)) # If the `reference_density_spacing` is set calculate the `ideal_neighbor_count` diff --git a/test/systems/edac_system.jl b/test/systems/edac_system.jl index 95549ed819..d4e89a52c4 100644 --- a/test/systems/edac_system.jl +++ b/test/systems/edac_system.jl @@ -141,6 +141,7 @@ │ ν₍EDAC₎: ………………………………………………………… ≈ 0.226 │ │ smoothing kernel: ………………………………… Val │ │ tansport velocity formulation: Nothing │ + │ average pressure reduction: ……… no │ │ acceleration: …………………………………………… [0.0, 0.0] │ │ surface tension: …………………………………… nothing │ │ surface normal method: …………………… nothing │ @@ -220,22 +221,28 @@ fluid = rectangular_patch(particle_spacing, (3, 3), seed=1) - system = EntropicallyDampedSPHSystem(fluid, smoothing_kernel, - transport_velocity=TransportVelocityAdami(0.0), - smoothing_length, 0.0) - semi = Semidiscretization(system) + transport_velocity = [nothing, TransportVelocityAdami(10000.0)] + names = ["No TVF", "TransportVelocityAdami"] + @testset "$(names[i])" for i in eachindex(transport_velocity) + system = EntropicallyDampedSPHSystem(fluid, smoothing_kernel, + transport_velocity=transport_velocity[i], + average_pressure_reduction=true, + smoothing_length, 0.0) + semi = Semidiscretization(system) - TrixiParticles.initialize_neighborhood_searches!(semi) + TrixiParticles.initialize_neighborhood_searches!(semi) - u_ode = vec(fluid.coordinates) - v_ode = vec(vcat(fluid.velocity, fluid.pressure')) + u_ode = vec(fluid.coordinates) + v_ode = vec(vcat(fluid.velocity, fluid.pressure')) - TrixiParticles.update_average_pressure!(system, system.transport_velocity, v_ode, - u_ode, semi) + TrixiParticles.update_average_pressure!(system, + system.average_pressure_reduction, + v_ode, u_ode, semi) - @test all(i -> system.cache.neighbor_counter[i] == nparticles(system), - nparticles(system)) - @test all(i -> isapprox(system.cache.pressure_average[i], -50.968532955185964), - nparticles(system)) + @test all(i -> system.cache.neighbor_counter[i] == nparticles(system), + nparticles(system)) + @test all(i -> isapprox(system.cache.pressure_average[i], -50.968532955185964), + nparticles(system)) + end end end From 8806ee25b3a62162e0dbee98b1bd705ddc5a83d2 Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Sat, 16 Aug 2025 19:57:35 +0300 Subject: [PATCH 2/3] Fix TGV validation --- .../validation_taylor_green_vortex_2d.jl | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/validation/taylor_green_vortex_2d/validation_taylor_green_vortex_2d.jl b/validation/taylor_green_vortex_2d/validation_taylor_green_vortex_2d.jl index 4543412306..849ddd38b6 100644 --- a/validation/taylor_green_vortex_2d/validation_taylor_green_vortex_2d.jl +++ b/validation/taylor_green_vortex_2d/validation_taylor_green_vortex_2d.jl @@ -17,6 +17,12 @@ reynolds_number = 100.0 density_calculators = [ContinuityDensity(), SummationDensity()] perturb_coordinates = [false, true] +# Define `average_pressure` for WCSPH, so that we can use the same code below for WCSPH +@inline function TrixiParticles.average_pressure(system::WeaklyCompressibleSPHSystem, + particle) + return zero(eltype(system)) +end + function compute_l1v_error(system, v_ode, u_ode, semi, t) v_analytical_avg = 0.0 v_avg = 0.0 From a8277e92273b12ebc01e89036f65e744839e95e3 Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Sun, 17 Aug 2025 14:09:55 +0300 Subject: [PATCH 3/3] Reformat --- src/schemes/fluid/entropically_damped_sph/system.jl | 2 +- test/systems/edac_system.jl | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/schemes/fluid/entropically_damped_sph/system.jl b/src/schemes/fluid/entropically_damped_sph/system.jl index d197e5ef02..19e50c441a 100644 --- a/src/schemes/fluid/entropically_damped_sph/system.jl +++ b/src/schemes/fluid/entropically_damped_sph/system.jl @@ -84,7 +84,7 @@ function EntropicallyDampedSPHSystem(initial_condition, smoothing_kernel, pressure_acceleration=inter_particle_averaged_pressure, density_calculator=SummationDensity(), transport_velocity=nothing, - average_pressure_reduction=!isnothing(transport_velocity), + average_pressure_reduction=(!isnothing(transport_velocity)), alpha=0.5, viscosity=nothing, acceleration=ntuple(_ -> 0.0, ndims(smoothing_kernel)), diff --git a/test/systems/edac_system.jl b/test/systems/edac_system.jl index d4e89a52c4..3d8cf9babb 100644 --- a/test/systems/edac_system.jl +++ b/test/systems/edac_system.jl @@ -240,9 +240,9 @@ v_ode, u_ode, semi) @test all(i -> system.cache.neighbor_counter[i] == nparticles(system), - nparticles(system)) + nparticles(system)) @test all(i -> isapprox(system.cache.pressure_average[i], -50.968532955185964), - nparticles(system)) + nparticles(system)) end end end