From 28f0c1937bdd2a58692994b7c6ad78b4db1bfddb Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Mon, 21 Jul 2025 15:38:36 +0200 Subject: [PATCH 01/21] Restructure TVF --- src/callbacks/update.jl | 4 - src/general/semidiscretization.jl | 13 +- src/general/system.jl | 4 + .../fluid/entropically_damped_sph/rhs.jl | 13 +- .../fluid/entropically_damped_sph/system.jl | 28 ++- src/schemes/fluid/fluid.jl | 23 -- src/schemes/fluid/transport_velocity.jl | 203 +++++++----------- .../fluid/weakly_compressible_sph/rhs.jl | 13 +- .../fluid/weakly_compressible_sph/system.jl | 28 ++- 9 files changed, 133 insertions(+), 196 deletions(-) diff --git a/src/callbacks/update.jl b/src/callbacks/update.jl index 7cc479d25b..44b24ed6ad 100644 --- a/src/callbacks/update.jl +++ b/src/callbacks/update.jl @@ -89,10 +89,6 @@ function (update_callback!::UpdateCallback)(integrator) update_particle_packing(system, v_ode, u_ode, semi, integrator) end - @trixi_timeit timer() "update TVF" foreach_system(semi) do system - update_transport_velocity!(system, v_ode, semi) - end - # Tell OrdinaryDiffEq that `u` has been modified u_modified!(integrator, true) diff --git a/src/general/semidiscretization.jl b/src/general/semidiscretization.jl index 3e033e447b..c4a9614fc4 100644 --- a/src/general/semidiscretization.jl +++ b/src/general/semidiscretization.jl @@ -464,8 +464,11 @@ function drift!(du_ode, v_ode, u_ode, semi, t) end @inline function add_velocity!(du, v, particle, system) + # This is zero unless a transport velocity is used + delta_v_ = delta_v(system, particle) + for i in 1:ndims(system) - du[i, particle] = v[i, particle] + du[i, particle] = v[i, particle] + delta_v_[i] end return du @@ -530,6 +533,14 @@ function update_systems_and_nhs(v_ode, u_ode, semi, t; update_from_callback=fals update_final!(system, v, u, v_ode, u_ode, semi, t; update_from_callback) end + + # Final update step for all remaining systems + foreach_system(semi) do system + v = wrap_v(v_ode, system, semi) + u = wrap_u(u_ode, system, semi) + + update_final2!(system, v, u, v_ode, u_ode, semi, t; update_from_callback) + end end function update_nhs!(semi, u_ode) diff --git a/src/general/system.jl b/src/general/system.jl index 69141a8ff0..869757f599 100644 --- a/src/general/system.jl +++ b/src/general/system.jl @@ -143,6 +143,10 @@ function update_final!(system, v, u, v_ode, u_ode, semi, t; update_from_callback return system end +function update_final2!(system, v, u, v_ode, u_ode, semi, t; update_from_callback=false) + return system +end + # Only for systems requiring a mandatory callback reset_callback_flag!(system) = system diff --git a/src/schemes/fluid/entropically_damped_sph/rhs.jl b/src/schemes/fluid/entropically_damped_sph/rhs.jl index 06e34272b7..e9e583ef46 100644 --- a/src/schemes/fluid/entropically_damped_sph/rhs.jl +++ b/src/schemes/fluid/entropically_damped_sph/rhs.jl @@ -50,10 +50,11 @@ function interact!(dv, v_particle_system, u_particle_system, sound_speed, m_a, m_b, rho_a, rho_b, grad_kernel) # Add convection term (only when using `TransportVelocityAdami`) - dv_convection = momentum_convection(particle_system, neighbor_system, - v_particle_system, v_neighbor_system, - rho_a, rho_b, m_a, m_b, - particle, neighbor, grad_kernel) + dv_convection = dv_transport_velocity(transport_velocity(particle_system), + particle_system, neighbor_system, + v_particle_system, v_neighbor_system, + rho_a, rho_b, m_a, m_b, + particle, neighbor, grad_kernel) dv_surface_tension = surface_tension_force(surface_tension_a, surface_tension_b, particle_system, neighbor_system, @@ -76,10 +77,6 @@ function interact!(dv, v_particle_system, u_particle_system, particle, neighbor, pos_diff, distance, sound_speed, m_a, m_b, p_a, p_b, rho_a, rho_b) - # Apply the transport velocity (only when using a transport velocity) - transport_velocity!(dv, particle_system, neighbor_system, particle, neighbor, - m_a, m_b, grad_kernel) - continuity_equation!(dv, density_calculator, v_diff, particle, m_b, rho_a, rho_b, particle_system, grad_kernel) end diff --git a/src/schemes/fluid/entropically_damped_sph/system.jl b/src/schemes/fluid/entropically_damped_sph/system.jl index 044104916b..af5594ea65 100644 --- a/src/schemes/fluid/entropically_damped_sph/system.jl +++ b/src/schemes/fluid/entropically_damped_sph/system.jl @@ -130,7 +130,7 @@ function EntropicallyDampedSPHSystem(initial_condition, smoothing_kernel, nu_edac = (alpha * smoothing_length * sound_speed) / 8 cache = (; create_cache_density(initial_condition, density_calculator)..., - create_cache_tvf(Val(:edac), initial_condition, transport_velocity)..., + create_cache_tvf_edac(initial_condition, transport_velocity)..., create_cache_surface_normal(surface_normal_method, ELTYPE, NDIMS, n_particles)..., create_cache_surface_tension(surface_tension, ELTYPE, NDIMS, @@ -208,14 +208,15 @@ function Base.show(io::IO, ::MIME"text/plain", system::EntropicallyDampedSPHSyst end end -create_cache_tvf(::Val{:edac}, initial_condition, ::Nothing) = (;) +create_cache_tvf_edac(initial_condition, ::Nothing) = (;) -function create_cache_tvf(::Val{:edac}, initial_condition, ::TransportVelocityAdami) +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)) update_callback_used = Ref(false) - return (; pressure_average, neighbor_counter, update_callback_used) + return (; delta_v, pressure_average, neighbor_counter, update_callback_used) end @inline function Base.eltype(::EntropicallyDampedSPHSystem{<:Any, ELTYPE}) where {ELTYPE} @@ -226,12 +227,12 @@ end return v_nvariables(system, system.density_calculator) end -@inline function v_nvariables(system::EntropicallyDampedSPHSystem, density_calculator) - return ndims(system) * factor_tvf(system) + 1 +@inline function v_nvariables(system::EntropicallyDampedSPHSystem, ::SummationDensity) + return ndims(system) + 1 end @inline function v_nvariables(system::EntropicallyDampedSPHSystem, ::ContinuityDensity) - return ndims(system) * factor_tvf(system) + 2 + return ndims(system) + 2 end system_correction(system::EntropicallyDampedSPHSystem) = system.correction @@ -311,10 +312,6 @@ function update_final!(system::EntropicallyDampedSPHSystem, v, u, v_ode, u_ode, compute_curvature!(system, surface_tension, v, u, v_ode, u_ode, semi, t) compute_stress_tensors!(system, surface_tension, v, u, v_ode, u_ode, semi, t) update_average_pressure!(system, system.transport_velocity, v_ode, u_ode, semi) - - # Check that TVF is only used together with `UpdateCallback` - check_tvf_configuration(system, system.transport_velocity, v, u, v_ode, u_ode, semi, t; - update_from_callback) end function update_average_pressure!(system, ::Nothing, v_ode, u_ode, semi) @@ -360,6 +357,15 @@ function update_average_pressure!(system, ::TransportVelocityAdami, v_ode, u_ode return system end +function update_final2!(system::EntropicallyDampedSPHSystem, v, u, v_ode, u_ode, semi, t; + update_from_callback=false) + if !update_from_callback + update_tvf!(system, transport_velocity(system), v, u, v_ode, u_ode, semi, t) + end + + return system +end + function write_v0!(v0, system::EntropicallyDampedSPHSystem, ::SummationDensity) # Note that `.=` is very slightly faster, but not GPU-compatible v0[end, :] = system.initial_condition.pressure diff --git a/src/schemes/fluid/fluid.jl b/src/schemes/fluid/fluid.jl index 308217170d..d3e6adf971 100644 --- a/src/schemes/fluid/fluid.jl +++ b/src/schemes/fluid/fluid.jl @@ -85,7 +85,6 @@ function write_v0!(v0, system::FluidSystem) copyto!(v0, indices, system.initial_condition.velocity, indices) write_v0!(v0, system, system.density_calculator) - write_v0!(v0, system, system.transport_velocity) return v0 end @@ -199,25 +198,3 @@ include("surface_tension.jl") include("surface_normal_sph.jl") include("weakly_compressible_sph/weakly_compressible_sph.jl") include("entropically_damped_sph/entropically_damped_sph.jl") - -@inline function add_velocity!(du, v, particle, - system::Union{EntropicallyDampedSPHSystem, - WeaklyCompressibleSPHSystem}) - add_velocity!(du, v, particle, system, system.transport_velocity) -end - -@inline function momentum_convection(system, neighbor_system, - v_particle_system, v_neighbor_system, rho_a, rho_b, - m_a, m_b, particle, neighbor, grad_kernel) - return zero(grad_kernel) -end - -@inline function momentum_convection(system, - neighbor_system::Union{EntropicallyDampedSPHSystem, - WeaklyCompressibleSPHSystem}, - v_particle_system, v_neighbor_system, rho_a, rho_b, - m_a, m_b, particle, neighbor, grad_kernel) - momentum_convection(system, neighbor_system, system.transport_velocity, - v_particle_system, v_neighbor_system, rho_a, rho_b, - m_a, m_b, particle, neighbor, grad_kernel) -end diff --git a/src/schemes/fluid/transport_velocity.jl b/src/schemes/fluid/transport_velocity.jl index 32ea57e7c9..284e7d876c 100644 --- a/src/schemes/fluid/transport_velocity.jl +++ b/src/schemes/fluid/transport_velocity.jl @@ -7,176 +7,117 @@ See [TVF](@ref transport_velocity_formulation) for more details of the method. # Arguments - `background_pressure`: Background pressure. Suggested is a background pressure which is on the order of the reference pressure. - -!!! note - There is no need for an artificial viscosity to suppress tensile instability when using `TransportVelocityAdami`. - Thus, it is highly recommended to use [`ViscosityAdami`](@ref) as viscosity model, - since [`ArtificialViscosityMonaghan`](@ref) leads to bad results. """ struct TransportVelocityAdami{T <: Real} background_pressure::T end -# Calculate `v_nvariables` appropriately -@inline factor_tvf(system::FluidSystem) = factor_tvf(system, system.transport_velocity) -@inline factor_tvf(system, ::Nothing) = 1 -@inline factor_tvf(system, ::TransportVelocityAdami) = 2 - -@inline update_transport_velocity!(system, v_ode, semi) = system +# No TVF for a system by default +@inline transport_velocity(system) = nothing -@inline function update_transport_velocity!(system::FluidSystem, v_ode, semi) - update_transport_velocity!(system, v_ode, semi, system.transport_velocity) +# `δ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) + return delta_v(system, transport_velocity(system), particle) end -@inline update_transport_velocity!(system, v_ode, semi, ::Nothing) = system - -@inline function update_transport_velocity!(system, v_ode, semi, ::TransportVelocityAdami) - v = wrap_v(v_ode, system, semi) - @threaded semi for particle in each_moving_particle(system) - for i in 1:ndims(system) - v[ndims(system) + i, particle] = v[i, particle] - end - end - - return system +@propagate_inbounds function delta_v(system, ::TransportVelocityAdami, particle) + return extract_svector(system.cache.delta_v, system, particle) end -function write_v0!(v0, system::FluidSystem, ::TransportVelocityAdami) - (; initial_condition) = system - - for particle in eachparticle(system) - # Write particle velocities - for dim in 1:ndims(system) - v0[ndims(system) + dim, particle] = initial_condition.velocity[dim, particle] - end - end - - return v0 +# Zero when no TVF is used +@inline function delta_v(system, transport_velocity, particle) + return zero(SVector{ndims(system), eltype(system)}) end -# Add momentum velocity. -@inline function add_velocity!(du, v, particle, system, ::Nothing) - for i in 1:ndims(system) - du[i, particle] = v[i, particle] - end - - return du +@inline function dv_transport_velocity(::Nothing, system, neighbor_system, + v_particle_system, v_neighbor_system, rho_a, rho_b, + m_a, m_b, particle, neighbor, grad_kernel) + return zero(grad_kernel) end -# Add advection velocity. -@inline function add_velocity!(du, v, particle, system, ::TransportVelocityAdami) - for i in 1:ndims(system) - du[i, particle] = v[ndims(system) + i, particle] - end +# @inline function dv_transport_velocity(::TransportVelocityAdami, system, neighbor_system, +# v_particle_system, v_neighbor_system, rho_a, rho_b, +# m_a, m_b, particle, neighbor, grad_kernel) +# v_a = current_velocity(v_particle_system, system, particle) +# delta_v_a = delta_v(system, particle) - return du -end +# v_b = current_velocity(v_neighbor_system, neighbor_system, neighbor) +# delta_v_b = delta_v(neighbor_system, neighbor) -@inline function advection_velocity(v, system, particle) - return SVector(ntuple(@inline(dim->v[ndims(system) + dim, particle]), ndims(system))) -end +# A_a = rho_a * v_a * delta_v_a' +# A_b = rho_b * v_b * delta_v_b' -@inline function momentum_convection(system, neighbor_system, ::Nothing, - v_particle_system, v_neighbor_system, rho_a, rho_b, - m_a, m_b, particle, neighbor, grad_kernel) - return zero(grad_kernel) -end +# return m_b * ((A_a / rho_a^2 + A_b / rho_b^2) / 2) * grad_kernel +# end -@inline function momentum_convection(system, neighbor_system, ::TransportVelocityAdami, - v_particle_system, v_neighbor_system, rho_a, rho_b, - m_a, m_b, particle, neighbor, grad_kernel) +@inline function dv_transport_velocity(::TransportVelocityAdami, system, neighbor_system, + v_particle_system, v_neighbor_system, rho_a, rho_b, + m_a, m_b, particle, neighbor, grad_kernel) volume_a = m_a / rho_a volume_b = m_b / rho_b volume_term = (volume_a^2 + volume_b^2) / m_a - momentum_velocity_a = current_velocity(v_particle_system, system, particle) - advection_velocity_a = advection_velocity(v_particle_system, system, particle) + v_a = current_velocity(v_particle_system, system, particle) + delta_v_a = delta_v(system, particle) - momentum_velocity_b = current_velocity(v_neighbor_system, neighbor_system, neighbor) - advection_velocity_b = advection_velocity(v_neighbor_system, neighbor_system, neighbor) + v_b = current_velocity(v_neighbor_system, neighbor_system, neighbor) + delta_v_b = delta_v(neighbor_system, neighbor) - A_a = rho_a * momentum_velocity_a * (advection_velocity_a - momentum_velocity_a)' - A_b = rho_b * momentum_velocity_b * (advection_velocity_b - momentum_velocity_b)' + A_a = rho_a * v_a * delta_v_a' + A_b = rho_b * v_b * delta_v_b' return volume_term * ((A_a + A_b) / 2) * grad_kernel end -@inline transport_velocity!(dv, system, neighbor_system, particle, neighbor, m_a, m_b, - grad_kernel) = dv - -@inline function transport_velocity!(dv, system::FluidSystem, neighbor_system, - particle, neighbor, m_a, m_b, grad_kernel) - transport_velocity!(dv, system, neighbor_system, system.transport_velocity, - particle, neighbor, m_a, m_b, grad_kernel) -end - -@inline function transport_velocity!(dv, system, neighbor_system, ::Nothing, - particle, neighbor, m_a, m_b, grad_kernel) - return dv +function update_tvf!(system, transport_velocity, v, u, v_ode, u_ode, semi, t) + return system end -@inline function transport_velocity!(dv, system, neighbor_system, ::TransportVelocityAdami, - particle, neighbor, m_a, m_b, grad_kernel) - (; transport_velocity) = system +function update_tvf!(system, transport_velocity::TransportVelocityAdami, v, u, v_ode, + u_ode, semi, t) + (; delta_v) = system.cache (; background_pressure) = transport_velocity - NDIMS = ndims(system) - - # The TVF is based on the assumption that the pressure gradient is only accurately - # computed when the particle distribution is isotropic. - # That means, the force contribution vanishes only if the particle distribution is - # isotropic AND the field being differentiated by the kernel gradient is spatially constant. - # So we must guarantee a constant field and therefore the reference density is used - # instead of the locally computed one. - # TODO: - # volume_a = particle_spacing(system, particle)^ndims(system) - # volume_b = particle_spacing(neighbor_system, neighbor)^ndims(neighbor_system) - volume_a = m_a / system.initial_condition.density[particle] - volume_b = m_b / neighbor_system.initial_condition.density[neighbor] - - volume_term = (volume_a^2 + volume_b^2) / m_a - - for dim in 1:NDIMS - dv[NDIMS + dim, particle] -= volume_term * background_pressure * grad_kernel[dim] - end - - return dv -end - -function reset_callback_flag!(system::FluidSystem) - reset_callback_flag!(system, system.transport_velocity) -end + set_zero!(delta_v) -reset_callback_flag!(system, ::Nothing) = system + foreach_system(semi) do neighbor_system + u_neighbor = wrap_u(u_ode, neighbor_system, semi) -function reset_callback_flag!(system::FluidSystem, ::TransportVelocityAdami) - system.cache.update_callback_used[] = false + system_coords = current_coordinates(u, system) + neighbor_coords = current_coordinates(u_neighbor, neighbor_system) - return system -end + foreach_point_neighbor(system, neighbor_system, system_coords, neighbor_coords, + semi; + points=each_moving_particle(system)) do particle, neighbor, + pos_diff, distance + m_a = hydrodynamic_mass(neighbor_system, neighbor) + m_b = hydrodynamic_mass(neighbor_system, neighbor) -function update_callback_used!(system::FluidSystem) - update_callback_used!(system, system.transport_velocity) -end + grad_kernel = smoothing_kernel_grad(system, pos_diff, distance, particle) -update_callback_used!(system, ::Nothing) = system + # The TVF is based on the assumption that the pressure gradient is only accurately + # computed when the particle distribution is isotropic. + # That means, the force contribution vanishes only if the particle distribution is + # isotropic AND the field being differentiated by the kernel gradient is spatially constant. + # So we must guarantee a constant field and therefore the reference density is used + # instead of the locally computed one. + # TODO: + # volume_a = particle_spacing(system, particle)^ndims(system) + # volume_b = particle_spacing(neighbor_system, neighbor)^ndims(neighbor_system) + volume_a = m_a / system.initial_condition.density[particle] + volume_b = m_b / neighbor_system.initial_condition.density[neighbor] -function update_callback_used!(system, transport_velocity) - system.cache.update_callback_used[] = true -end + volume_term = (volume_a^2 + volume_b^2) / m_a -function check_tvf_configuration(system::FluidSystem, ::Nothing, - v, u, v_ode, u_ode, semi, t; update_from_callback=false) - return system -end + delta_v_ = -volume_term * background_pressure * grad_kernel -function check_tvf_configuration(system::FluidSystem, ::TransportVelocityAdami, - v, u, v_ode, u_ode, semi, t; update_from_callback=false) - # The `UpdateCallback` will set `system.cache.update_callback_used[]` to `true` - # in the initialization. However, other callbacks might update the system first, hence `update_from_callback`. - if !update_from_callback && !(system.cache.update_callback_used[]) - throw(ArgumentError("`UpdateCallback` is required when using `TransportVelocityAdami`")) - end + # Write into the buffer + for i in eachindex(delta_v_) + @inbounds delta_v[i, particle] += delta_v_[i] + end + end + end - return system + return system end diff --git a/src/schemes/fluid/weakly_compressible_sph/rhs.jl b/src/schemes/fluid/weakly_compressible_sph/rhs.jl index a1987a8c0b..a42b0270a6 100644 --- a/src/schemes/fluid/weakly_compressible_sph/rhs.jl +++ b/src/schemes/fluid/weakly_compressible_sph/rhs.jl @@ -70,10 +70,11 @@ function interact!(dv, v_particle_system, u_particle_system, grad_kernel) # Add convection term (only when using `TransportVelocityAdami`) - dv_convection = momentum_convection(particle_system, neighbor_system, - v_particle_system, v_neighbor_system, - rho_a, rho_b, m_a, m_b, - particle, neighbor, grad_kernel) + dv_convection = dv_transport_velocity(transport_velocity(particle_system), + particle_system, neighbor_system, + v_particle_system, v_neighbor_system, + rho_a, rho_b, m_a, m_b, + particle, neighbor, grad_kernel) dv_surface_tension = surface_tension_correction * surface_tension_force(surface_tension_a, surface_tension_b, @@ -93,10 +94,6 @@ function interact!(dv, v_particle_system, u_particle_system, # debug_array[i, particle] += dv_pressure[i] end - # Apply the transport velocity (only when using a transport velocity) - transport_velocity!(dv, particle_system, neighbor_system, particle, neighbor, - m_a, m_b, grad_kernel) - # TODO If variable smoothing_length is used, this should use the neighbor smoothing length # Propagate `@inbounds` to the continuity equation, which accesses particle data @inbounds continuity_equation!(dv, density_calculator, v_particle_system, diff --git a/src/schemes/fluid/weakly_compressible_sph/system.jl b/src/schemes/fluid/weakly_compressible_sph/system.jl index 690803960f..543e3f3d45 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(Val(:wcsph), initial_condition, transport_velocity)..., + create_cache_tvf_wcsph(initial_condition, transport_velocity)..., color=Int(color_value)) # If the `reference_density_spacing` is set calculate the `ideal_neighbor_count` @@ -221,12 +221,13 @@ function Base.show(io::IO, ::MIME"text/plain", system::WeaklyCompressibleSPHSyst end end -create_cache_tvf(::Val{:wcsph}, initial_condition, ::Nothing) = (;) +create_cache_tvf_wcsph(initial_condition, ::Nothing) = (;) -function create_cache_tvf(::Val{:wcsph}, initial_condition, ::TransportVelocityAdami) +function create_cache_tvf_wcsph(initial_condition, ::TransportVelocityAdami) + delta_v = zeros(eltype(initial_condition), ndims(initial_condition), nparticles(initial_condition)) update_callback_used = Ref(false) - return (; update_callback_used) + return (; delta_v, update_callback_used) end @inline function Base.eltype(::WeaklyCompressibleSPHSystem{<:Any, ELTYPE}) where {ELTYPE} @@ -237,12 +238,12 @@ end return v_nvariables(system, system.density_calculator) end -@inline function v_nvariables(system::WeaklyCompressibleSPHSystem, density_calculator) - return ndims(system) * factor_tvf(system) +@inline function v_nvariables(system::WeaklyCompressibleSPHSystem, ::SummationDensity) + return ndims(system) end @inline function v_nvariables(system::WeaklyCompressibleSPHSystem, ::ContinuityDensity) - return ndims(system) * factor_tvf(system) + 1 + return ndims(system) + 1 end system_correction(system::WeaklyCompressibleSPHSystem) = system.correction @@ -286,6 +287,8 @@ end @inline system_sound_speed(system::WeaklyCompressibleSPHSystem) = system.state_equation.sound_speed +@inline transport_velocity(system::WeaklyCompressibleSPHSystem) = system.transport_velocity + function update_quantities!(system::WeaklyCompressibleSPHSystem, v, u, v_ode, u_ode, semi, t) (; density_calculator, density_diffusion, correction) = system @@ -321,10 +324,15 @@ function update_final!(system::WeaklyCompressibleSPHSystem, 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) +end - # Check that TVF is only used together with `UpdateCallback` - check_tvf_configuration(system, system.transport_velocity, v, u, v_ode, u_ode, semi, t; - update_from_callback) +function update_final2!(system::WeaklyCompressibleSPHSystem, v, u, v_ode, u_ode, semi, t; + update_from_callback=false) + if !update_from_callback + update_tvf!(system, transport_velocity(system), v, u, v_ode, u_ode, semi, t) + end + + return system end function kernel_correct_density!(system::WeaklyCompressibleSPHSystem, v, u, v_ode, u_ode, From f42b0faa623bbe535e4973cccf000c91c94170fe Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Thu, 31 Jul 2025 16:05:20 +0200 Subject: [PATCH 02/21] Remove dv term from newer paper --- src/schemes/fluid/transport_velocity.jl | 15 --------------- 1 file changed, 15 deletions(-) diff --git a/src/schemes/fluid/transport_velocity.jl b/src/schemes/fluid/transport_velocity.jl index 284e7d876c..74a0f2f356 100644 --- a/src/schemes/fluid/transport_velocity.jl +++ b/src/schemes/fluid/transport_velocity.jl @@ -36,21 +36,6 @@ end return zero(grad_kernel) end -# @inline function dv_transport_velocity(::TransportVelocityAdami, system, neighbor_system, -# v_particle_system, v_neighbor_system, rho_a, rho_b, -# m_a, m_b, particle, neighbor, grad_kernel) -# v_a = current_velocity(v_particle_system, system, particle) -# delta_v_a = delta_v(system, particle) - -# v_b = current_velocity(v_neighbor_system, neighbor_system, neighbor) -# delta_v_b = delta_v(neighbor_system, neighbor) - -# A_a = rho_a * v_a * delta_v_a' -# A_b = rho_b * v_b * delta_v_b' - -# return m_b * ((A_a / rho_a^2 + A_b / rho_b^2) / 2) * grad_kernel -# end - @inline function dv_transport_velocity(::TransportVelocityAdami, system, neighbor_system, v_particle_system, v_neighbor_system, rho_a, rho_b, m_a, m_b, particle, neighbor, grad_kernel) From cdb0770bff8abf5f2fee1bc65c841b3321e33d48 Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Thu, 31 Jul 2025 18:18:31 +0200 Subject: [PATCH 03/21] Eliminate time step from spatial discretization --- src/schemes/fluid/transport_velocity.jl | 97 +++++++++++++------ .../fluid/weakly_compressible_sph/rhs.jl | 8 +- 2 files changed, 72 insertions(+), 33 deletions(-) diff --git a/src/schemes/fluid/transport_velocity.jl b/src/schemes/fluid/transport_velocity.jl index 74a0f2f356..c76e18772a 100644 --- a/src/schemes/fluid/transport_velocity.jl +++ b/src/schemes/fluid/transport_velocity.jl @@ -31,19 +31,17 @@ end end @inline function dv_transport_velocity(::Nothing, system, neighbor_system, - v_particle_system, v_neighbor_system, rho_a, rho_b, - m_a, m_b, particle, neighbor, grad_kernel) + particle, neighbor, v_system, v_neighbor_system, + m_a, m_b, rho_a, rho_b, pos_diff, distance, + grad_kernel, correction) return zero(grad_kernel) end @inline function dv_transport_velocity(::TransportVelocityAdami, system, neighbor_system, - v_particle_system, v_neighbor_system, rho_a, rho_b, - m_a, m_b, particle, neighbor, grad_kernel) - volume_a = m_a / rho_a - volume_b = m_b / rho_b - volume_term = (volume_a^2 + volume_b^2) / m_a - - v_a = current_velocity(v_particle_system, system, particle) + particle, neighbor, v_system, v_neighbor_system, + m_a, m_b, rho_a, rho_b, pos_diff, distance, + grad_kernel, correction) + v_a = current_velocity(v_system, system, particle) delta_v_a = delta_v(system, particle) v_b = current_velocity(v_neighbor_system, neighbor_system, neighbor) @@ -52,7 +50,16 @@ end A_a = rho_a * v_a * delta_v_a' A_b = rho_b * v_b * delta_v_b' - return volume_term * ((A_a + A_b) / 2) * grad_kernel + # The following term depends on the pressure acceleration formulation. + # See the large comment below. In the original paper (Adami et al., 2013), this is + # (V_a^2 + V_b^2) / m_a * ((A_a + A_b) / 2) * ∇W_ab. + # With the most common pressure acceleration formulation, this is + # m_b * (A_a + A_b) / (ρ_a * ρ_b) * ∇W_ab. + # In order to obtain this, we pass `p_a = A_a` and `p_b = A_b` to the + # `pressure_acceleration` function. + return pressure_acceleration(system, neighbor_system, particle, neighbor, + m_a, m_b, A_a, A_b, rho_a, rho_b, pos_diff, + distance, grad_kernel, correction) end function update_tvf!(system, transport_velocity, v, u, v_ode, u_ode, semi, t) @@ -61,9 +68,12 @@ end function update_tvf!(system, transport_velocity::TransportVelocityAdami, v, u, v_ode, u_ode, semi, t) - (; delta_v) = system.cache + (; cache, correction) = system + (; delta_v) = cache (; background_pressure) = transport_velocity + sound_speed = system_sound_speed(system) + set_zero!(delta_v) foreach_system(semi) do neighbor_system @@ -76,33 +86,60 @@ function update_tvf!(system, transport_velocity::TransportVelocityAdami, v, u, v semi; points=each_moving_particle(system)) do particle, neighbor, pos_diff, distance - m_a = hydrodynamic_mass(neighbor_system, neighbor) - m_b = hydrodynamic_mass(neighbor_system, neighbor) + m_a = @inbounds hydrodynamic_mass(neighbor_system, neighbor) + m_b = @inbounds hydrodynamic_mass(neighbor_system, neighbor) - grad_kernel = smoothing_kernel_grad(system, pos_diff, distance, particle) + rho_a = @inbounds current_density(system, particle) + rho_b = @inbounds current_density(neighbor_system, neighbor) - # The TVF is based on the assumption that the pressure gradient is only accurately - # computed when the particle distribution is isotropic. - # That means, the force contribution vanishes only if the particle distribution is - # isotropic AND the field being differentiated by the kernel gradient is spatially constant. - # So we must guarantee a constant field and therefore the reference density is used - # instead of the locally computed one. - # TODO: - # volume_a = particle_spacing(system, particle)^ndims(system) - # volume_b = particle_spacing(neighbor_system, neighbor)^ndims(neighbor_system) - volume_a = m_a / system.initial_condition.density[particle] - volume_b = m_b / neighbor_system.initial_condition.density[neighbor] + h = smoothing_length(system, particle) - volume_term = (volume_a^2 + volume_b^2) / m_a + grad_kernel = smoothing_kernel_grad(system, pos_diff, distance, particle) - delta_v_ = -volume_term * background_pressure * grad_kernel + # In the original paper (Adami et al., 2013), the transport velocity is applied + # as follows: + # v_{1/2} = v_0 + Δt/2 * a, + # where a is the regular SPH acceleration term (pressure, viscosity, etc.). + # r_1 = r_0 + Δt * (v_{1/2} + δv_{1/2}), + # where v_{1/2} + δv_{1/2} is the transport velocity. + # We will call δv_{1/2} the shifting velocity, which is given by + # δv = -Δt * p_0 / m_a * \sum_b (V_a^2 + V_b^2) * ∇W_ab, + # where p_0 is the background pressure, V_a = m_a / ρ_a, V_b = m_b / ρ_b. + # This term depends on the pressure acceleration formulation. + # In Zhang et al. (2017), the pressure acceleration term + # m_b * (p_a / ρ_a^2 + p_b / ρ_b^2) * ∇W_ab + # is used. They consequently changed the shifting velocity to + # δv = -Δt * p_0 * m_b * (1 / ρ_a^2 + 1 / ρ_b^2) * ∇W_ab. + # We therefore use the function `pressure_acceleration` to compute the + # shifting velocity according to the used pressure acceleration formulation. + # In most cases, this will be + # δv = -Δt * p_0 * m_b * (1 + 1) / (ρ_a * ρ_b^2) * ∇W_ab. + # + # In these papers, the shifting velocity is scaled by the time step Δt. + # We generally want the spatial discretization to be independent of the time step. + # Scaling the shifting velocity by the time step would lead to less shifting + # when very small time steps are used for testing/debugging purposes. + # This is especially problematic in TrixiParticles.jl, as the time step can vary + # significantly between different time integration methods (low vs high order). + # In order to eliminate the time step from the shifting velocity, we apply the + # CFL condition used in Adami et al. (2013): + # Δt <= 0.25 * h / c, + # where h is the smoothing length and c is the sound speed. + # Applying this equation as equality yields the shifting velocity + # δv = -p_0 / 4 * h / c * m_b * (1 + 1) / (ρ_a * ρ_b^2) * ∇W_ab. + # The last part is achieved by passing `p_a = 1` and `p_b = 1` to the + # `pressure_acceleration` function. + delta_v_ = -background_pressure / 4 * h / sound_speed * + pressure_acceleration(system, neighbor_system, particle, neighbor, + m_a, m_b, 1, 1, rho_a, rho_b, pos_diff, + distance, grad_kernel, correction) # Write into the buffer for i in eachindex(delta_v_) @inbounds delta_v[i, particle] += delta_v_[i] end - end - end + end + end - return system + return system end diff --git a/src/schemes/fluid/weakly_compressible_sph/rhs.jl b/src/schemes/fluid/weakly_compressible_sph/rhs.jl index a42b0270a6..03e94d2694 100644 --- a/src/schemes/fluid/weakly_compressible_sph/rhs.jl +++ b/src/schemes/fluid/weakly_compressible_sph/rhs.jl @@ -6,7 +6,8 @@ function interact!(dv, v_particle_system, u_particle_system, v_neighbor_system, u_neighbor_system, particle_system::WeaklyCompressibleSPHSystem, neighbor_system, semi) (; density_calculator, state_equation, correction) = particle_system - (; sound_speed) = state_equation + + sound_speed = system_sound_speed(particle_system) surface_tension_a = surface_tension_model(particle_system) surface_tension_b = surface_tension_model(neighbor_system) @@ -72,9 +73,10 @@ function interact!(dv, v_particle_system, u_particle_system, # Add convection term (only when using `TransportVelocityAdami`) dv_convection = dv_transport_velocity(transport_velocity(particle_system), particle_system, neighbor_system, + particle, neighbor, v_particle_system, v_neighbor_system, - rho_a, rho_b, m_a, m_b, - particle, neighbor, grad_kernel) + m_a, m_b, rho_a, rho_b, pos_diff, distance, + grad_kernel, correction) dv_surface_tension = surface_tension_correction * surface_tension_force(surface_tension_a, surface_tension_b, From 4ea73958de54dcb3c5aaffd20120f62687ca9750 Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Thu, 31 Jul 2025 18:19:46 +0200 Subject: [PATCH 04/21] Rename updates --- src/general/semidiscretization.jl | 6 +++--- src/general/system.jl | 5 +++-- .../open_boundary/method_of_characteristics.jl | 3 ++- src/schemes/boundary/open_boundary/mirroring.jl | 2 +- src/schemes/boundary/open_boundary/system.jl | 6 +++--- src/schemes/boundary/system.jl | 6 +++--- src/schemes/fluid/entropically_damped_sph/system.jl | 10 +--------- src/schemes/fluid/weakly_compressible_sph/system.jl | 10 +--------- src/schemes/solid/total_lagrangian_sph/system.jl | 4 ++-- 9 files changed, 19 insertions(+), 33 deletions(-) diff --git a/src/general/semidiscretization.jl b/src/general/semidiscretization.jl index c4a9614fc4..6074c6a2fb 100644 --- a/src/general/semidiscretization.jl +++ b/src/general/semidiscretization.jl @@ -526,12 +526,12 @@ function update_systems_and_nhs(v_ode, u_ode, semi, t; update_from_callback=fals update_pressure!(system, v, u, v_ode, u_ode, semi, t) end - # Final update step for all remaining systems foreach_system(semi) do system v = wrap_v(v_ode, system, semi) u = wrap_u(u_ode, system, semi) - update_final!(system, v, u, v_ode, u_ode, semi, t; update_from_callback) + update_boundary_interpolation!(system, v, u, v_ode, u_ode, semi, t; + update_from_callback) end # Final update step for all remaining systems @@ -539,7 +539,7 @@ function update_systems_and_nhs(v_ode, u_ode, semi, t; update_from_callback=fals v = wrap_v(v_ode, system, semi) u = wrap_u(u_ode, system, semi) - update_final2!(system, v, u, v_ode, u_ode, semi, t; update_from_callback) + update_final!(system, v, u, v_ode, u_ode, semi, t; update_from_callback) end end diff --git a/src/general/system.jl b/src/general/system.jl index 869757f599..ed1bba4092 100644 --- a/src/general/system.jl +++ b/src/general/system.jl @@ -139,11 +139,12 @@ function update_pressure!(system, v, u, v_ode, u_ode, semi, t) return system end -function update_final!(system, v, u, v_ode, u_ode, semi, t; update_from_callback=false) +function update_boundary_interpolation!(system, v, u, v_ode, u_ode, semi, t; + update_from_callback=false) return system end -function update_final2!(system, v, u, v_ode, u_ode, semi, t; update_from_callback=false) +function update_final!(system, v, u, v_ode, u_ode, semi, t; update_from_callback=false) return system end diff --git a/src/schemes/boundary/open_boundary/method_of_characteristics.jl b/src/schemes/boundary/open_boundary/method_of_characteristics.jl index 8dacc446e7..67c4f47d6c 100644 --- a/src/schemes/boundary/open_boundary/method_of_characteristics.jl +++ b/src/schemes/boundary/open_boundary/method_of_characteristics.jl @@ -75,7 +75,8 @@ end end # Called from semidiscretization -function update_final!(system, ::BoundaryModelLastiwka, v, u, v_ode, u_ode, semi, t) +function update_boundary_model!(system, ::BoundaryModelLastiwka, v, u, v_ode, u_ode, + semi, t) @trixi_timeit timer() "evaluate characteristics" begin evaluate_characteristics!(system, v, u, v_ode, u_ode, semi, t) end diff --git a/src/schemes/boundary/open_boundary/mirroring.jl b/src/schemes/boundary/open_boundary/mirroring.jl index 2eb7b286de..0ea726da95 100644 --- a/src/schemes/boundary/open_boundary/mirroring.jl +++ b/src/schemes/boundary/open_boundary/mirroring.jl @@ -24,7 +24,7 @@ function update_boundary_quantities!(system, ::BoundaryModelTafuni, v, u, v_ode, end end -update_final!(system, ::BoundaryModelTafuni, v, u, v_ode, u_ode, semi, t) = system +update_boundary_model!(system, ::BoundaryModelTafuni, v, u, v_ode, u_ode, semi, t) = system function extrapolate_values!(system, v_open_boundary, v_fluid, u_open_boundary, u_fluid, semi, t; prescribed_density=false, diff --git a/src/schemes/boundary/open_boundary/system.jl b/src/schemes/boundary/open_boundary/system.jl index 14034fedc2..4b722ad90c 100644 --- a/src/schemes/boundary/open_boundary/system.jl +++ b/src/schemes/boundary/open_boundary/system.jl @@ -255,13 +255,13 @@ end return system.pressure end -function update_final!(system::OpenBoundarySPHSystem, v, u, v_ode, u_ode, semi, t; - update_from_callback=false) +function update_boundary_interpolation!(system::OpenBoundarySPHSystem, v, u, v_ode, u_ode, + semi, t; update_from_callback=false) if !update_from_callback && !(system.update_callback_used[]) throw(ArgumentError("`UpdateCallback` is required when using `OpenBoundarySPHSystem`")) end - update_final!(system, system.boundary_model, v, u, v_ode, u_ode, semi, t) + update_boundary_model!(system, system.boundary_model, v, u, v_ode, u_ode, semi, t) end # This function is called by the `UpdateCallback`, as the integrator array might be modified diff --git a/src/schemes/boundary/system.jl b/src/schemes/boundary/system.jl index 4b01a403ec..fb1b2e01dd 100644 --- a/src/schemes/boundary/system.jl +++ b/src/schemes/boundary/system.jl @@ -358,9 +358,9 @@ function update_quantities!(system::BoundarySPHSystem, v, u, v_ode, u_ode, semi, end # This update depends on the computed quantities of the fluid system and therefore -# has to be in `update_final!` after `update_quantities!`. -function update_final!(system::BoundarySPHSystem, v, u, v_ode, u_ode, semi, t; - update_from_callback=false) +# has to be in `update_boundary_interpolation!` after `update_quantities!`. +function update_boundary_interpolation!(system::BoundarySPHSystem, v, u, v_ode, u_ode, semi, + t; update_from_callback=false) (; boundary_model) = system # Note that `update_pressure!(::BoundarySPHSystem, ...)` is empty diff --git a/src/schemes/fluid/entropically_damped_sph/system.jl b/src/schemes/fluid/entropically_damped_sph/system.jl index af5594ea65..258342bf32 100644 --- a/src/schemes/fluid/entropically_damped_sph/system.jl +++ b/src/schemes/fluid/entropically_damped_sph/system.jl @@ -312,6 +312,7 @@ function update_final!(system::EntropicallyDampedSPHSystem, v, u, v_ode, u_ode, compute_curvature!(system, surface_tension, v, u, v_ode, u_ode, semi, t) compute_stress_tensors!(system, surface_tension, v, u, v_ode, u_ode, semi, t) update_average_pressure!(system, system.transport_velocity, 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) @@ -357,15 +358,6 @@ function update_average_pressure!(system, ::TransportVelocityAdami, v_ode, u_ode return system end -function update_final2!(system::EntropicallyDampedSPHSystem, v, u, v_ode, u_ode, semi, t; - update_from_callback=false) - if !update_from_callback - update_tvf!(system, transport_velocity(system), v, u, v_ode, u_ode, semi, t) - end - - return system -end - function write_v0!(v0, system::EntropicallyDampedSPHSystem, ::SummationDensity) # Note that `.=` is very slightly faster, but not GPU-compatible v0[end, :] = system.initial_condition.pressure diff --git a/src/schemes/fluid/weakly_compressible_sph/system.jl b/src/schemes/fluid/weakly_compressible_sph/system.jl index 543e3f3d45..eea22ab9d3 100644 --- a/src/schemes/fluid/weakly_compressible_sph/system.jl +++ b/src/schemes/fluid/weakly_compressible_sph/system.jl @@ -324,15 +324,7 @@ function update_final!(system::WeaklyCompressibleSPHSystem, 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) -end - -function update_final2!(system::WeaklyCompressibleSPHSystem, v, u, v_ode, u_ode, semi, t; - update_from_callback=false) - if !update_from_callback - update_tvf!(system, transport_velocity(system), v, u, v_ode, u_ode, semi, t) - end - - return system + update_tvf!(system, transport_velocity(system), v, u, v_ode, u_ode, semi, t) end function kernel_correct_density!(system::WeaklyCompressibleSPHSystem, v, u, v_ode, u_ode, diff --git a/src/schemes/solid/total_lagrangian_sph/system.jl b/src/schemes/solid/total_lagrangian_sph/system.jl index acf269d7ba..7440d70763 100644 --- a/src/schemes/solid/total_lagrangian_sph/system.jl +++ b/src/schemes/solid/total_lagrangian_sph/system.jl @@ -261,8 +261,8 @@ function update_quantities!(system::TotalLagrangianSPHSystem, v, u, v_ode, u_ode return system end -function update_final!(system::TotalLagrangianSPHSystem, v, u, v_ode, u_ode, semi, t; - update_from_callback=false) +function update_boundary_interpolation!(system::TotalLagrangianSPHSystem, v, u, + v_ode, u_ode, semi, t; update_from_callback=false) (; boundary_model) = system # Only update boundary model From 4aaf014065650ca4ed74ff318bdb2b449cb027bb Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Thu, 31 Jul 2025 18:39:08 +0200 Subject: [PATCH 05/21] Fix EDAC --- src/schemes/fluid/entropically_damped_sph/rhs.jl | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/schemes/fluid/entropically_damped_sph/rhs.jl b/src/schemes/fluid/entropically_damped_sph/rhs.jl index e9e583ef46..548b1ff4c4 100644 --- a/src/schemes/fluid/entropically_damped_sph/rhs.jl +++ b/src/schemes/fluid/entropically_damped_sph/rhs.jl @@ -52,9 +52,10 @@ function interact!(dv, v_particle_system, u_particle_system, # Add convection term (only when using `TransportVelocityAdami`) dv_convection = dv_transport_velocity(transport_velocity(particle_system), particle_system, neighbor_system, + particle, neighbor, v_particle_system, v_neighbor_system, - rho_a, rho_b, m_a, m_b, - particle, neighbor, grad_kernel) + m_a, m_b, rho_a, rho_b, pos_diff, distance, + grad_kernel, correction) dv_surface_tension = surface_tension_force(surface_tension_a, surface_tension_b, particle_system, neighbor_system, From 16f38d2dbb26d484abae7080db8a8518509b50d3 Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Thu, 31 Jul 2025 18:40:51 +0200 Subject: [PATCH 06/21] Fix `current_density` --- src/schemes/fluid/transport_velocity.jl | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/schemes/fluid/transport_velocity.jl b/src/schemes/fluid/transport_velocity.jl index c76e18772a..33199c1961 100644 --- a/src/schemes/fluid/transport_velocity.jl +++ b/src/schemes/fluid/transport_velocity.jl @@ -77,6 +77,7 @@ function update_tvf!(system, transport_velocity::TransportVelocityAdami, v, u, v set_zero!(delta_v) foreach_system(semi) do neighbor_system + v_neighbor = wrap_v(v_ode, neighbor_system, semi) u_neighbor = wrap_u(u_ode, neighbor_system, semi) system_coords = current_coordinates(u, system) @@ -89,8 +90,8 @@ function update_tvf!(system, transport_velocity::TransportVelocityAdami, v, u, v m_a = @inbounds hydrodynamic_mass(neighbor_system, neighbor) m_b = @inbounds hydrodynamic_mass(neighbor_system, neighbor) - rho_a = @inbounds current_density(system, particle) - rho_b = @inbounds current_density(neighbor_system, neighbor) + rho_a = @inbounds current_density(v, system, particle) + rho_b = @inbounds current_density(v_neighbor, neighbor_system, neighbor) h = smoothing_length(system, particle) From f9456d681cf4574f5763326f19f6d986c11c1cfe Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Thu, 31 Jul 2025 18:44:35 +0200 Subject: [PATCH 07/21] Fix sign of background pressure --- src/schemes/fluid/entropically_damped_sph/system.jl | 3 ++- src/schemes/fluid/transport_velocity.jl | 2 +- src/schemes/fluid/weakly_compressible_sph/system.jl | 3 ++- 3 files changed, 5 insertions(+), 3 deletions(-) diff --git a/src/schemes/fluid/entropically_damped_sph/system.jl b/src/schemes/fluid/entropically_damped_sph/system.jl index 258342bf32..58b39b0c63 100644 --- a/src/schemes/fluid/entropically_damped_sph/system.jl +++ b/src/schemes/fluid/entropically_damped_sph/system.jl @@ -211,7 +211,8 @@ 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)) + 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)) update_callback_used = Ref(false) diff --git a/src/schemes/fluid/transport_velocity.jl b/src/schemes/fluid/transport_velocity.jl index c76e18772a..11ca5537e9 100644 --- a/src/schemes/fluid/transport_velocity.jl +++ b/src/schemes/fluid/transport_velocity.jl @@ -129,7 +129,7 @@ function update_tvf!(system, transport_velocity::TransportVelocityAdami, v, u, v # δv = -p_0 / 4 * h / c * m_b * (1 + 1) / (ρ_a * ρ_b^2) * ∇W_ab. # The last part is achieved by passing `p_a = 1` and `p_b = 1` to the # `pressure_acceleration` function. - delta_v_ = -background_pressure / 4 * h / sound_speed * + delta_v_ = background_pressure / 4 * h / sound_speed * pressure_acceleration(system, neighbor_system, particle, neighbor, m_a, m_b, 1, 1, rho_a, rho_b, pos_diff, distance, grad_kernel, correction) diff --git a/src/schemes/fluid/weakly_compressible_sph/system.jl b/src/schemes/fluid/weakly_compressible_sph/system.jl index eea22ab9d3..79990be2c4 100644 --- a/src/schemes/fluid/weakly_compressible_sph/system.jl +++ b/src/schemes/fluid/weakly_compressible_sph/system.jl @@ -224,7 +224,8 @@ end create_cache_tvf_wcsph(initial_condition, ::Nothing) = (;) function create_cache_tvf_wcsph(initial_condition, ::TransportVelocityAdami) - delta_v = zeros(eltype(initial_condition), ndims(initial_condition), nparticles(initial_condition)) + delta_v = zeros(eltype(initial_condition), ndims(initial_condition), + nparticles(initial_condition)) update_callback_used = Ref(false) return (; delta_v, update_callback_used) From 3157178574a3427a57d3d2f4db61acd63e24e39c Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Thu, 31 Jul 2025 18:48:20 +0200 Subject: [PATCH 08/21] Add reference --- src/schemes/fluid/transport_velocity.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/schemes/fluid/transport_velocity.jl b/src/schemes/fluid/transport_velocity.jl index a178c81144..e0d0eb463f 100644 --- a/src/schemes/fluid/transport_velocity.jl +++ b/src/schemes/fluid/transport_velocity.jl @@ -1,7 +1,8 @@ """ TransportVelocityAdami(background_pressure::Real) -Transport Velocity Formulation (TVF) to suppress pairing and tensile instability. +Transport Velocity Formulation (TVF) by [Adami et al. (2013)](@cite Adami2013) +to suppress pairing and tensile instability. See [TVF](@ref transport_velocity_formulation) for more details of the method. # Arguments From 5725d90540e81fed998cdcfa281dae099c1faf39 Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Fri, 1 Aug 2025 14:21:13 +0200 Subject: [PATCH 09/21] Implement suggestions --- src/general/semidiscretization.jl | 2 ++ src/general/system.jl | 2 +- .../fluid/entropically_damped_sph/rhs.jl | 18 +++++++++--------- src/schemes/fluid/transport_velocity.jl | 2 +- .../fluid/weakly_compressible_sph/rhs.jl | 17 ++++++++--------- 5 files changed, 21 insertions(+), 20 deletions(-) diff --git a/src/general/semidiscretization.jl b/src/general/semidiscretization.jl index 6074c6a2fb..8dbcab6e34 100644 --- a/src/general/semidiscretization.jl +++ b/src/general/semidiscretization.jl @@ -526,6 +526,8 @@ function update_systems_and_nhs(v_ode, u_ode, semi, t; update_from_callback=fals update_pressure!(system, v, u, v_ode, u_ode, semi, t) end + # This update depends on the computed quantities of the fluid system and therefore + # needs to be after `update_quantities!`. foreach_system(semi) do system v = wrap_v(v_ode, system, semi) u = wrap_u(u_ode, system, semi) diff --git a/src/general/system.jl b/src/general/system.jl index ed1bba4092..118ded4c94 100644 --- a/src/general/system.jl +++ b/src/general/system.jl @@ -126,7 +126,7 @@ end system_correction(system), system, particle) end -# System update orders. This can be dispatched if needed. +# System updates do nothing by default, but can be dispatched if needed function update_positions!(system, v, u, v_ode, u_ode, semi, t) return system end diff --git a/src/schemes/fluid/entropically_damped_sph/rhs.jl b/src/schemes/fluid/entropically_damped_sph/rhs.jl index 548b1ff4c4..7a0d1ded58 100644 --- a/src/schemes/fluid/entropically_damped_sph/rhs.jl +++ b/src/schemes/fluid/entropically_damped_sph/rhs.jl @@ -50,12 +50,12 @@ function interact!(dv, v_particle_system, u_particle_system, sound_speed, m_a, m_b, rho_a, rho_b, grad_kernel) # Add convection term (only when using `TransportVelocityAdami`) - dv_convection = dv_transport_velocity(transport_velocity(particle_system), - particle_system, neighbor_system, - particle, neighbor, - v_particle_system, v_neighbor_system, - m_a, m_b, rho_a, rho_b, pos_diff, distance, - grad_kernel, correction) + dv_tvf = dv_transport_velocity(transport_velocity(particle_system), + particle_system, neighbor_system, + particle, neighbor, + v_particle_system, v_neighbor_system, + m_a, m_b, rho_a, rho_b, pos_diff, distance, + grad_kernel, correction) dv_surface_tension = surface_tension_force(surface_tension_a, surface_tension_b, particle_system, neighbor_system, @@ -66,9 +66,9 @@ function interact!(dv, v_particle_system, u_particle_system, particle, neighbor, pos_diff, distance) for i in 1:ndims(particle_system) - dv[i, - particle] += dv_pressure[i] + dv_viscosity_[i] + dv_convection[i] + - dv_surface_tension[i] + dv_adhesion[i] + @inbounds dv[i, + particle] += dv_pressure[i] + dv_viscosity_[i] + dv_tvf[i] + + dv_surface_tension[i] + dv_adhesion[i] end v_diff = current_velocity(v_particle_system, particle_system, particle) - diff --git a/src/schemes/fluid/transport_velocity.jl b/src/schemes/fluid/transport_velocity.jl index e0d0eb463f..18dd49da92 100644 --- a/src/schemes/fluid/transport_velocity.jl +++ b/src/schemes/fluid/transport_velocity.jl @@ -128,7 +128,7 @@ function update_tvf!(system, transport_velocity::TransportVelocityAdami, v, u, v # Δt <= 0.25 * h / c, # where h is the smoothing length and c is the sound speed. # Applying this equation as equality yields the shifting velocity - # δv = -p_0 / 4 * h / c * m_b * (1 + 1) / (ρ_a * ρ_b^2) * ∇W_ab. + # δv = -p_0 / 4 * h / c * m_b * (1 + 1) / (ρ_a * ρ_b) * ∇W_ab. # The last part is achieved by passing `p_a = 1` and `p_b = 1` to the # `pressure_acceleration` function. delta_v_ = background_pressure / 4 * h / sound_speed * diff --git a/src/schemes/fluid/weakly_compressible_sph/rhs.jl b/src/schemes/fluid/weakly_compressible_sph/rhs.jl index 03e94d2694..ecc6975d94 100644 --- a/src/schemes/fluid/weakly_compressible_sph/rhs.jl +++ b/src/schemes/fluid/weakly_compressible_sph/rhs.jl @@ -71,12 +71,12 @@ function interact!(dv, v_particle_system, u_particle_system, grad_kernel) # Add convection term (only when using `TransportVelocityAdami`) - dv_convection = dv_transport_velocity(transport_velocity(particle_system), - particle_system, neighbor_system, - particle, neighbor, - v_particle_system, v_neighbor_system, - m_a, m_b, rho_a, rho_b, pos_diff, distance, - grad_kernel, correction) + dv_tvf = dv_transport_velocity(transport_velocity(particle_system), + particle_system, neighbor_system, + particle, neighbor, + v_particle_system, v_neighbor_system, + m_a, m_b, rho_a, rho_b, pos_diff, distance, + grad_kernel, correction) dv_surface_tension = surface_tension_correction * surface_tension_force(surface_tension_a, surface_tension_b, @@ -89,9 +89,8 @@ function interact!(dv, v_particle_system, u_particle_system, for i in 1:ndims(particle_system) @inbounds dv[i, - particle] += dv_pressure[i] + dv_viscosity_[i] + - dv_convection[i] + dv_surface_tension[i] + - dv_adhesion[i] + particle] += dv_pressure[i] + dv_viscosity_[i] + dv_tvf[i] + + dv_surface_tension[i] + dv_adhesion[i] # Debug example # debug_array[i, particle] += dv_pressure[i] end From 42fc6d67625985698e840b7696628b6382c272cc Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Fri, 1 Aug 2025 15:21:15 +0200 Subject: [PATCH 10/21] Remove `update_from_callback` --- src/callbacks/density_reinit.jl | 2 +- src/callbacks/particle_shifting.jl | 3 +-- src/callbacks/post_process.jl | 3 +-- src/callbacks/update.jl | 6 ++--- src/general/semidiscretization.jl | 27 ++++++++++++++----- src/general/system.jl | 19 ++++++------- src/io/write_vtk.jl | 3 +-- src/preprocessing/particle_packing/system.jl | 18 ++++--------- src/schemes/boundary/open_boundary/system.jl | 15 +++++------ src/schemes/boundary/system.jl | 3 +-- .../fluid/entropically_damped_sph/system.jl | 6 ++--- .../fluid/weakly_compressible_sph/system.jl | 6 ++--- .../solid/total_lagrangian_sph/system.jl | 2 +- 13 files changed, 53 insertions(+), 60 deletions(-) diff --git a/src/callbacks/density_reinit.jl b/src/callbacks/density_reinit.jl index 50f54e2860..ded47149f9 100644 --- a/src/callbacks/density_reinit.jl +++ b/src/callbacks/density_reinit.jl @@ -69,7 +69,7 @@ function initialize_reinit_cb!(cb::DensityReinitializationCallback, u, t, integr # Update systems to compute quantities like density and pressure. semi = integrator.p v_ode, u_ode = u.x - update_systems_and_nhs(v_ode, u_ode, semi, t; update_from_callback=true) + update_systems_and_nhs(v_ode, u_ode, semi, t) # Apply the callback. cb(integrator) diff --git a/src/callbacks/particle_shifting.jl b/src/callbacks/particle_shifting.jl index 5ff735fcb5..921370f3b1 100644 --- a/src/callbacks/particle_shifting.jl +++ b/src/callbacks/particle_shifting.jl @@ -39,8 +39,7 @@ function particle_shifting!(integrator) # still have the values from the last stage of the previous step if not updated here. @trixi_timeit timer() "update systems and nhs" begin # Don't create sub-timers here to avoid cluttering the timer output - @notimeit timer() update_systems_and_nhs(v_ode, u_ode, semi, t; - update_from_callback=true) + @notimeit timer() update_systems_and_nhs(v_ode, u_ode, semi, t) end @trixi_timeit timer() "particle shifting" foreach_system(semi) do system diff --git a/src/callbacks/post_process.jl b/src/callbacks/post_process.jl index 8764013166..36e16ce7f0 100644 --- a/src/callbacks/post_process.jl +++ b/src/callbacks/post_process.jl @@ -240,8 +240,7 @@ function (pp::PostprocessCallback)(integrator) # still have the values from the last stage of the previous step if not updated here. @trixi_timeit timer() "update systems and nhs" begin # Don't create sub-timers here to avoid cluttering the timer output - @notimeit timer() update_systems_and_nhs(v_ode, u_ode, semi, t; - update_from_callback=true) + @notimeit timer() update_systems_and_nhs(v_ode, u_ode, semi, t) end foreach_system(semi) do system diff --git a/src/callbacks/update.jl b/src/callbacks/update.jl index 44b24ed6ad..c3d5c59298 100644 --- a/src/callbacks/update.jl +++ b/src/callbacks/update.jl @@ -57,7 +57,7 @@ function initial_update!(cb::UpdateCallback, u, t, integrator) # Tell systems that `UpdateCallback` is used foreach_system(semi) do system - update_callback_used!(system) + set_callback_flag!(system, true) end return cb(integrator) @@ -78,7 +78,7 @@ function (update_callback!::UpdateCallback)(integrator) # Update quantities that are stored in the systems. These quantities (e.g. pressure) # still have the values from the last stage of the previous step if not updated here. - update_systems_and_nhs(v_ode, u_ode, semi, t; update_from_callback=true) + update_systems_and_nhs(v_ode, u_ode, semi, t) # Update open boundaries first, since particles might be activated or deactivated @trixi_timeit timer() "update open boundary" foreach_system(semi) do system @@ -137,5 +137,3 @@ function Base.show(io::IO, ::MIME"text/plain", summary_box(io, "UpdateCallback", setup) end end - -update_callback_used!(system) = system diff --git a/src/general/semidiscretization.jl b/src/general/semidiscretization.jl index 8dbcab6e34..7c22d541e1 100644 --- a/src/general/semidiscretization.jl +++ b/src/general/semidiscretization.jl @@ -338,7 +338,7 @@ function semidiscretize(semi, tspan; reset_threads=true) # Initialize this system initialize!(system, semi_new) - # Only for systems requiring a mandatory callback + # Only for systems requiring the use of the `UpdateCallback` reset_callback_flag!(system) end @@ -373,7 +373,7 @@ function restart_with!(semi, sol; reset_threads=true) restart_with!(system, v, u) - # Only for systems requiring a mandatory callback + # Only for systems requiring the use of the `UpdateCallback` reset_callback_flag!(system) end @@ -478,6 +478,9 @@ end function kick!(dv_ode, v_ode, u_ode, semi, t) @trixi_timeit timer() "kick!" begin + # Check that the `UpdateCallback` is used if required + check_update_callback(semi) + @trixi_timeit timer() "reset ∂v/∂t" set_zero!(dv_ode) @trixi_timeit timer() "update systems and nhs" update_systems_and_nhs(v_ode, u_ode, @@ -493,8 +496,9 @@ function kick!(dv_ode, v_ode, u_ode, semi, t) return dv_ode end -# Update the systems and neighborhood searches (NHS) for a simulation before calling `interact!` to compute forces -function update_systems_and_nhs(v_ode, u_ode, semi, t; update_from_callback=false) +# Update the systems and neighborhood searches (NHS) for a simulation +# before calling `interact!` to compute forces. +function update_systems_and_nhs(v_ode, u_ode, semi, t) # First update step before updating the NHS # (for example for writing the current coordinates in the solid system) foreach_system(semi) do system @@ -532,8 +536,7 @@ function update_systems_and_nhs(v_ode, u_ode, semi, t; update_from_callback=fals v = wrap_v(v_ode, system, semi) u = wrap_u(u_ode, system, semi) - update_boundary_interpolation!(system, v, u, v_ode, u_ode, semi, t; - update_from_callback) + update_boundary_interpolation!(system, v, u, v_ode, u_ode, semi, t) end # Final update step for all remaining systems @@ -541,7 +544,7 @@ function update_systems_and_nhs(v_ode, u_ode, semi, t; update_from_callback=fals v = wrap_v(v_ode, system, semi) u = wrap_u(u_ode, system, semi) - update_final!(system, v, u, v_ode, u_ode, semi, t; update_from_callback) + update_final!(system, v, u, v_ode, u_ode, semi, t) end end @@ -871,6 +874,16 @@ function update!(neighborhood_search, x, y, semi; points_moving=(true, false), parallelization_backend=semi.parallelization_backend) end +function check_update_callback(semi) + foreach_system(semi) do system + # This check will be optimized away if the system does not require the callback + if requires_update_callback(system) && !update_callback_used(system) + system_name = system |> typeof |> nameof + throw(ArgumentError("`UpdateCallback` is required for `$system_name`")) + end + end +end + function check_configuration(systems, nhs::Union{Nothing, PointNeighbors.AbstractNeighborhoodSearch}) foreach_system(systems) do system diff --git a/src/general/system.jl b/src/general/system.jl index 118ded4c94..c757617719 100644 --- a/src/general/system.jl +++ b/src/general/system.jl @@ -139,26 +139,27 @@ function update_pressure!(system, v, u, v_ode, u_ode, semi, t) return system end -function update_boundary_interpolation!(system, v, u, v_ode, u_ode, semi, t; - update_from_callback=false) +function update_boundary_interpolation!(system, v, u, v_ode, u_ode, semi, t) return system end -function update_final!(system, v, u, v_ode, u_ode, semi, t; update_from_callback=false) +function update_final!(system, v, u, v_ode, u_ode, semi, t) return system end -# Only for systems requiring a mandatory callback -reset_callback_flag!(system) = system +# Only for systems requiring the use of the `UpdateCallback` +@inline requires_update_callback(system) = false +@inline update_callback_used(system) = false +@inline set_callback_flag!(system, value) = system -initial_smoothing_length(system) = smoothing_length(system, nothing) +@inline initial_smoothing_length(system) = smoothing_length(system, nothing) -function smoothing_length(system, particle) +@inline function smoothing_length(system, particle) return system.smoothing_length end -system_smoothing_kernel(system) = system.smoothing_kernel -system_correction(system) = nothing +@inline system_smoothing_kernel(system) = system.smoothing_kernel +@inline system_correction(system) = nothing @inline particle_spacing(system, particle) = system.initial_condition.particle_spacing diff --git a/src/io/write_vtk.jl b/src/io/write_vtk.jl index a2b40f31e1..b93e25f324 100644 --- a/src/io/write_vtk.jl +++ b/src/io/write_vtk.jl @@ -58,8 +58,7 @@ function trixi2vtk(vu_ode, semi, t; iter=nothing, output_directory="out", prefix # still have the values from the last stage of the previous step if not updated here. @trixi_timeit timer() "update systems" begin # Don't create sub-timers here to avoid cluttering the timer output - @notimeit timer() update_systems_and_nhs(v_ode, u_ode, semi, t; - update_from_callback=true) + @notimeit timer() update_systems_and_nhs(v_ode, u_ode, semi, t) end filenames = system_names(systems) diff --git a/src/preprocessing/particle_packing/system.jl b/src/preprocessing/particle_packing/system.jl index 1318f63544..0e149178ea 100644 --- a/src/preprocessing/particle_packing/system.jl +++ b/src/preprocessing/particle_packing/system.jl @@ -215,14 +215,15 @@ end return ndims(system) end -function reset_callback_flag!(system::ParticlePackingSystem) - system.update_callback_used[] = false +@inline requires_update_callback(system::ParticlePackingSystem) = true +@inline update_callback_used(system::ParticlePackingSystem) = system.update_callback_used[] + +function set_callback_flag!(system::ParticlePackingSystem, value) + system.update_callback_used[] = value return system end -update_callback_used!(system::ParticlePackingSystem) = system.update_callback_used[] = true - function write2vtk!(vtk, v, u, t, system::ParticlePackingSystem; write_meta_data=true) vtk["velocity"] = [advection_velocity(v, system, particle) for particle in active_particles(system)] @@ -289,15 +290,6 @@ function update_position!(u, system::ParticlePackingSystem, semi) return u end -function update_final!(system::ParticlePackingSystem, v, u, v_ode, u_ode, semi, t; - update_from_callback=false) - if !update_from_callback && !(system.update_callback_used[]) - throw(ArgumentError("`UpdateCallback` is required when using `ParticlePackingSystem`")) - end - - return system -end - # Skip for systems without `SignedDistanceField` constrain_particles_onto_surface!(u, system::ParticlePackingSystem{Nothing}, semi) = u diff --git a/src/schemes/boundary/open_boundary/system.jl b/src/schemes/boundary/open_boundary/system.jl index 60155662c9..e820bd7aef 100644 --- a/src/schemes/boundary/open_boundary/system.jl +++ b/src/schemes/boundary/open_boundary/system.jl @@ -226,14 +226,15 @@ end return ELTYPE end -function reset_callback_flag!(system::OpenBoundarySPHSystem) - system.update_callback_used[] = false +@inline requires_update_callback(system::OpenBoundarySPHSystem) = true +@inline update_callback_used(system::OpenBoundarySPHSystem) = system.update_callback_used[] + +function set_callback_flag!(system::OpenBoundarySPHSystem, value) + system.update_callback_used[] = value return system end -update_callback_used!(system::OpenBoundarySPHSystem) = system.update_callback_used[] = true - function corresponding_fluid_system(system::OpenBoundarySPHSystem, semi) return system.fluid_system end @@ -253,11 +254,7 @@ end end function update_boundary_interpolation!(system::OpenBoundarySPHSystem, v, u, v_ode, u_ode, - semi, t; update_from_callback=false) - if !update_from_callback && !(system.update_callback_used[]) - throw(ArgumentError("`UpdateCallback` is required when using `OpenBoundarySPHSystem`")) - end - + semi) update_boundary_model!(system, system.boundary_model, v, u, v_ode, u_ode, semi, t) end diff --git a/src/schemes/boundary/system.jl b/src/schemes/boundary/system.jl index fb1b2e01dd..84e2611a8b 100644 --- a/src/schemes/boundary/system.jl +++ b/src/schemes/boundary/system.jl @@ -359,8 +359,7 @@ end # This update depends on the computed quantities of the fluid system and therefore # has to be in `update_boundary_interpolation!` after `update_quantities!`. -function update_boundary_interpolation!(system::BoundarySPHSystem, v, u, v_ode, u_ode, semi, - t; update_from_callback=false) +function update_boundary_interpolation!(system::BoundarySPHSystem, v, u, v_ode, u_ode, semi) (; boundary_model) = system # Note that `update_pressure!(::BoundarySPHSystem, ...)` is empty diff --git a/src/schemes/fluid/entropically_damped_sph/system.jl b/src/schemes/fluid/entropically_damped_sph/system.jl index 58b39b0c63..49e670cba4 100644 --- a/src/schemes/fluid/entropically_damped_sph/system.jl +++ b/src/schemes/fluid/entropically_damped_sph/system.jl @@ -215,9 +215,8 @@ function create_cache_tvf_edac(initial_condition, ::TransportVelocityAdami) nparticles(initial_condition)) pressure_average = copy(initial_condition.pressure) neighbor_counter = Vector{Int}(undef, nparticles(initial_condition)) - update_callback_used = Ref(false) - return (; delta_v, pressure_average, neighbor_counter, update_callback_used) + return (; delta_v, pressure_average, neighbor_counter) end @inline function Base.eltype(::EntropicallyDampedSPHSystem{<:Any, ELTYPE}) where {ELTYPE} @@ -305,8 +304,7 @@ function update_pressure!(system::EntropicallyDampedSPHSystem, v, u, v_ode, u_od compute_surface_delta_function!(system, system.surface_tension, semi) end -function update_final!(system::EntropicallyDampedSPHSystem, v, u, v_ode, u_ode, semi, t; - update_from_callback=false) +function update_final!(system::EntropicallyDampedSPHSystem, v, u, v_ode, u_ode, semi, t) (; surface_tension) = system # Surface normal of neighbor and boundary needs to have been calculated already diff --git a/src/schemes/fluid/weakly_compressible_sph/system.jl b/src/schemes/fluid/weakly_compressible_sph/system.jl index 79990be2c4..c8d72fd5a0 100644 --- a/src/schemes/fluid/weakly_compressible_sph/system.jl +++ b/src/schemes/fluid/weakly_compressible_sph/system.jl @@ -226,9 +226,8 @@ create_cache_tvf_wcsph(initial_condition, ::Nothing) = (;) function create_cache_tvf_wcsph(initial_condition, ::TransportVelocityAdami) delta_v = zeros(eltype(initial_condition), ndims(initial_condition), nparticles(initial_condition)) - update_callback_used = Ref(false) - return (; delta_v, update_callback_used) + return (; delta_v) end @inline function Base.eltype(::WeaklyCompressibleSPHSystem{<:Any, ELTYPE}) where {ELTYPE} @@ -318,8 +317,7 @@ function update_pressure!(system::WeaklyCompressibleSPHSystem, v, u, v_ode, u_od return system end -function update_final!(system::WeaklyCompressibleSPHSystem, v, u, v_ode, u_ode, semi, t; - update_from_callback=false) +function update_final!(system::WeaklyCompressibleSPHSystem, v, u, v_ode, u_ode, semi, t) (; surface_tension) = system # Surface normal of neighbor and boundary needs to have been calculated already diff --git a/src/schemes/solid/total_lagrangian_sph/system.jl b/src/schemes/solid/total_lagrangian_sph/system.jl index 7440d70763..b465ca2d76 100644 --- a/src/schemes/solid/total_lagrangian_sph/system.jl +++ b/src/schemes/solid/total_lagrangian_sph/system.jl @@ -262,7 +262,7 @@ function update_quantities!(system::TotalLagrangianSPHSystem, v, u, v_ode, u_ode end function update_boundary_interpolation!(system::TotalLagrangianSPHSystem, v, u, - v_ode, u_ode, semi, t; update_from_callback=false) + v_ode, u_ode, semi, t) (; boundary_model) = system # Only update boundary model From 23c63d796b89ea9376642ffaeeb7e5e6922e9a7c Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Fri, 1 Aug 2025 15:24:04 +0200 Subject: [PATCH 11/21] Implement suggestions --- .../fluid/entropically_damped_sph/system.jl | 13 +----------- src/schemes/fluid/transport_velocity.jl | 20 +++++++++++++++++++ .../fluid/weakly_compressible_sph/system.jl | 9 --------- 3 files changed, 21 insertions(+), 21 deletions(-) diff --git a/src/schemes/fluid/entropically_damped_sph/system.jl b/src/schemes/fluid/entropically_damped_sph/system.jl index 49e670cba4..45c5e8074d 100644 --- a/src/schemes/fluid/entropically_damped_sph/system.jl +++ b/src/schemes/fluid/entropically_damped_sph/system.jl @@ -208,17 +208,6 @@ function Base.show(io::IO, ::MIME"text/plain", system::EntropicallyDampedSPHSyst end 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 - @inline function Base.eltype(::EntropicallyDampedSPHSystem{<:Any, ELTYPE}) where {ELTYPE} return ELTYPE end @@ -264,7 +253,7 @@ end @inline average_pressure(system, particle) = zero(eltype(system)) @inline function average_pressure(system::EntropicallyDampedSPHSystem, particle) - average_pressure(system, system.transport_velocity, particle) + average_pressure(system, transport_velocity(system), particle) end @inline function average_pressure(system, ::TransportVelocityAdami, particle) diff --git a/src/schemes/fluid/transport_velocity.jl b/src/schemes/fluid/transport_velocity.jl index 18dd49da92..d7a978e356 100644 --- a/src/schemes/fluid/transport_velocity.jl +++ b/src/schemes/fluid/transport_velocity.jl @@ -16,6 +16,26 @@ end # No TVF for a system by default @inline transport_velocity(system) = nothing +create_cache_tvf_wcsph(initial_condition, ::Nothing) = (;) + +function create_cache_tvf_wcsph(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 c8d72fd5a0..e8de016119 100644 --- a/src/schemes/fluid/weakly_compressible_sph/system.jl +++ b/src/schemes/fluid/weakly_compressible_sph/system.jl @@ -221,15 +221,6 @@ function Base.show(io::IO, ::MIME"text/plain", system::WeaklyCompressibleSPHSyst end end -create_cache_tvf_wcsph(initial_condition, ::Nothing) = (;) - -function create_cache_tvf_wcsph(initial_condition, ::TransportVelocityAdami) - delta_v = zeros(eltype(initial_condition), ndims(initial_condition), - nparticles(initial_condition)) - - return (; delta_v) -end - @inline function Base.eltype(::WeaklyCompressibleSPHSystem{<:Any, ELTYPE}) where {ELTYPE} return ELTYPE end From 49417b537a38ff596595be3eaf9b6ded8a9d2233 Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Fri, 1 Aug 2025 15:38:17 +0200 Subject: [PATCH 12/21] Fix tests --- src/general/semidiscretization.jl | 4 +-- .../fluid/entropically_damped_sph/system.jl | 7 ++-- test/schemes/fluid/fluid.jl | 1 - test/schemes/fluid/transport_velocity.jl | 35 ------------------- test/systems/edac_system.jl | 2 +- 5 files changed, 7 insertions(+), 42 deletions(-) delete mode 100644 test/schemes/fluid/transport_velocity.jl diff --git a/src/general/semidiscretization.jl b/src/general/semidiscretization.jl index 7c22d541e1..70f5510b65 100644 --- a/src/general/semidiscretization.jl +++ b/src/general/semidiscretization.jl @@ -339,7 +339,7 @@ function semidiscretize(semi, tspan; reset_threads=true) initialize!(system, semi_new) # Only for systems requiring the use of the `UpdateCallback` - reset_callback_flag!(system) + set_callback_flag!(system, false) end return DynamicalODEProblem(kick!, drift!, v0_ode, u0_ode, tspan, semi_new) @@ -374,7 +374,7 @@ function restart_with!(semi, sol; reset_threads=true) restart_with!(system, v, u) # Only for systems requiring the use of the `UpdateCallback` - reset_callback_flag!(system) + set_callback_flag!(system, false) end return semi diff --git a/src/schemes/fluid/entropically_damped_sph/system.jl b/src/schemes/fluid/entropically_damped_sph/system.jl index 45c5e8074d..f9875e684a 100644 --- a/src/schemes/fluid/entropically_damped_sph/system.jl +++ b/src/schemes/fluid/entropically_damped_sph/system.jl @@ -250,6 +250,8 @@ end @inline system_sound_speed(system::EntropicallyDampedSPHSystem) = system.sound_speed +@inline transport_velocity(system::EntropicallyDampedSPHSystem) = system.transport_velocity + @inline average_pressure(system, particle) = zero(eltype(system)) @inline function average_pressure(system::EntropicallyDampedSPHSystem, particle) @@ -299,7 +301,7 @@ 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, system.transport_velocity, v_ode, u_ode, semi) + update_average_pressure!(system, transport_velocity(system), v_ode, u_ode, semi) update_tvf!(system, transport_velocity(system), v, u, v_ode, u_ode, semi, t) end @@ -333,8 +335,7 @@ function update_average_pressure!(system, ::TransportVelocityAdami, v_ode, u_ode points=each_moving_particle(system)) do particle, neighbor, pos_diff, distance pressure_average[particle] += current_pressure(v_neighbor_system, - neighbor_system, - neighbor) + neighbor_system, neighbor) neighbor_counter[particle] += 1 end end diff --git a/test/schemes/fluid/fluid.jl b/test/schemes/fluid/fluid.jl index 570c3bc865..664325a1ad 100644 --- a/test/schemes/fluid/fluid.jl +++ b/test/schemes/fluid/fluid.jl @@ -2,6 +2,5 @@ include("weakly_compressible_sph/weakly_compressible_sph.jl") include("rhs.jl") include("pressure_acceleration.jl") include("surface_normal_sph.jl") -include("transport_velocity.jl") include("surface_tension.jl") include("viscosity.jl") diff --git a/test/schemes/fluid/transport_velocity.jl b/test/schemes/fluid/transport_velocity.jl deleted file mode 100644 index bfa71e1252..0000000000 --- a/test/schemes/fluid/transport_velocity.jl +++ /dev/null @@ -1,35 +0,0 @@ -@trixi_testset "Transport Velocity Formulation" begin - particle_spacing = 0.1 - smoothing_kernel = SchoenbergCubicSplineKernel{2}() - smoothing_length = 1.2particle_spacing - - fluid = rectangular_patch(particle_spacing, (3, 3)) - - v0_tvf = zeros(5, nparticles(fluid)) - - system_tvf = EntropicallyDampedSPHSystem(fluid, smoothing_kernel, - transport_velocity=TransportVelocityAdami(0.0), - smoothing_length, 0.0) - system = EntropicallyDampedSPHSystem(fluid, smoothing_kernel, smoothing_length, 0.0) - - @testset "Number of Variables" begin - @test TrixiParticles.v_nvariables(system_tvf) == 5 - @test TrixiParticles.v_nvariables(system) == 3 - end - - @testset "write_v0!" begin - TrixiParticles.write_v0!(v0_tvf, system_tvf) - - @test isapprox(vcat(fluid.velocity, fluid.velocity, fluid.pressure'), v0_tvf) - end - - @testset "Update" begin - semi = Semidiscretization(system_tvf) - fill!(v0_tvf, 1.5) - v0_tvf[1:2, :] .= 2.5 - - TrixiParticles.update_transport_velocity!(system_tvf, vec(v0_tvf), semi) - - @test isapprox(fill(2.5, (4, nparticles(system_tvf))), v0_tvf[1:4, :]) - end -end diff --git a/test/systems/edac_system.jl b/test/systems/edac_system.jl index 424680abf3..95549ed819 100644 --- a/test/systems/edac_system.jl +++ b/test/systems/edac_system.jl @@ -228,7 +228,7 @@ TrixiParticles.initialize_neighborhood_searches!(semi) u_ode = vec(fluid.coordinates) - v_ode = vec(vcat(fluid.velocity, fluid.velocity, fluid.pressure')) + v_ode = vec(vcat(fluid.velocity, fluid.pressure')) TrixiParticles.update_average_pressure!(system, system.transport_velocity, v_ode, u_ode, semi) From 11978d172c03adc8257f665a05220a3d318cb3c9 Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Fri, 1 Aug 2025 18:10:50 +0200 Subject: [PATCH 13/21] Fix background pressure factor --- src/schemes/fluid/transport_velocity.jl | 16 +++++++++------- 1 file changed, 9 insertions(+), 7 deletions(-) diff --git a/src/schemes/fluid/transport_velocity.jl b/src/schemes/fluid/transport_velocity.jl index d7a978e356..bdfa0af8f1 100644 --- a/src/schemes/fluid/transport_velocity.jl +++ b/src/schemes/fluid/transport_velocity.jl @@ -122,20 +122,22 @@ function update_tvf!(system, transport_velocity::TransportVelocityAdami, v, u, v # as follows: # v_{1/2} = v_0 + Δt/2 * a, # where a is the regular SPH acceleration term (pressure, viscosity, etc.). - # r_1 = r_0 + Δt * (v_{1/2} + δv_{1/2}), - # where v_{1/2} + δv_{1/2} is the transport velocity. + # r_1 = r_0 + Δt * (v_{1/2}, + # where ̃v_{1/2} = v_{1/2} + Δt/2 * p_0 / m_a * \sum_b[ (V_a^2 + V_b^2) * ∇W_ab ] + # is the transport velocity. + # We call δv_{1/2} = ̃v_{1/2} - v_{1/2} the shifting velocity. # We will call δv_{1/2} the shifting velocity, which is given by - # δv = -Δt * p_0 / m_a * \sum_b (V_a^2 + V_b^2) * ∇W_ab, + # δv = -Δt/2 * p_0 / m_a * \sum_b[ (V_a^2 + V_b^2) * ∇W_ab ], # where p_0 is the background pressure, V_a = m_a / ρ_a, V_b = m_b / ρ_b. # This term depends on the pressure acceleration formulation. # In Zhang et al. (2017), the pressure acceleration term # m_b * (p_a / ρ_a^2 + p_b / ρ_b^2) * ∇W_ab # is used. They consequently changed the shifting velocity to - # δv = -Δt * p_0 * m_b * (1 / ρ_a^2 + 1 / ρ_b^2) * ∇W_ab. + # δv = -Δt/2 * p_0 * \sum_b[ m_b * (1 / ρ_a^2 + 1 / ρ_b^2) * ∇W_ab ]. # We therefore use the function `pressure_acceleration` to compute the # shifting velocity according to the used pressure acceleration formulation. # In most cases, this will be - # δv = -Δt * p_0 * m_b * (1 + 1) / (ρ_a * ρ_b^2) * ∇W_ab. + # δv = -Δt/2 * p_0 * \sum_b[ m_b * (1 + 1) / (ρ_a * ρ_b) * ∇W_ab ]. # # In these papers, the shifting velocity is scaled by the time step Δt. # We generally want the spatial discretization to be independent of the time step. @@ -148,10 +150,10 @@ function update_tvf!(system, transport_velocity::TransportVelocityAdami, v, u, v # Δt <= 0.25 * h / c, # where h is the smoothing length and c is the sound speed. # Applying this equation as equality yields the shifting velocity - # δv = -p_0 / 4 * h / c * m_b * (1 + 1) / (ρ_a * ρ_b) * ∇W_ab. + # δv = -p_0 / 8 * h / c * \sum_b[ m_b * (1 + 1) / (ρ_a * ρ_b) * ∇W_ab ]. # The last part is achieved by passing `p_a = 1` and `p_b = 1` to the # `pressure_acceleration` function. - delta_v_ = background_pressure / 4 * h / sound_speed * + delta_v_ = background_pressure / 8 * h / sound_speed * pressure_acceleration(system, neighbor_system, particle, neighbor, m_a, m_b, 1, 1, rho_a, rho_b, pos_diff, distance, grad_kernel, correction) From a1344f9c81ea98b1350c9cda8d47f4d4439350e1 Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Fri, 1 Aug 2025 18:42:10 +0200 Subject: [PATCH 14/21] Improve comments --- src/general/semidiscretization.jl | 2 +- src/schemes/boundary/dummy_particles/dummy_particles.jl | 7 +++---- src/schemes/boundary/system.jl | 3 ++- src/schemes/fluid/weakly_compressible_sph/rhs.jl | 2 +- src/schemes/fluid/weakly_compressible_sph/system.jl | 8 +++++--- 5 files changed, 12 insertions(+), 10 deletions(-) diff --git a/src/general/semidiscretization.jl b/src/general/semidiscretization.jl index 70f5510b65..abe0d087f1 100644 --- a/src/general/semidiscretization.jl +++ b/src/general/semidiscretization.jl @@ -468,7 +468,7 @@ end delta_v_ = delta_v(system, particle) for i in 1:ndims(system) - du[i, particle] = v[i, particle] + delta_v_[i] + @inbounds du[i, particle] = v[i, particle] + delta_v_[i] end return du diff --git a/src/schemes/boundary/dummy_particles/dummy_particles.jl b/src/schemes/boundary/dummy_particles/dummy_particles.jl index 1d30a5274a..f0051bdb3a 100644 --- a/src/schemes/boundary/dummy_particles/dummy_particles.jl +++ b/src/schemes/boundary/dummy_particles/dummy_particles.jl @@ -311,17 +311,16 @@ end system, v, u, v_ode, u_ode, semi) (; correction, density_calculator) = boundary_model - compute_correction_values!(system, correction, u, v_ode, u_ode, semi) + compute_pressure!(boundary_model, density_calculator, system, v, u, v_ode, u_ode, semi) + # These are only computed when using corrections + compute_correction_values!(system, correction, u, v_ode, u_ode, semi) compute_gradient_correction_matrix!(correction, boundary_model, system, u, v_ode, u_ode, semi) - # `kernel_correct_density!` only performed for `SummationDensity` kernel_correct_density!(boundary_model, v, u, v_ode, u_ode, semi, correction, density_calculator) - compute_pressure!(boundary_model, density_calculator, system, v, u, v_ode, u_ode, semi) - return boundary_model end diff --git a/src/schemes/boundary/system.jl b/src/schemes/boundary/system.jl index 84e2611a8b..77607118f3 100644 --- a/src/schemes/boundary/system.jl +++ b/src/schemes/boundary/system.jl @@ -362,7 +362,8 @@ end function update_boundary_interpolation!(system::BoundarySPHSystem, v, u, v_ode, u_ode, semi) (; boundary_model) = system - # Note that `update_pressure!(::BoundarySPHSystem, ...)` is empty + # Note that `update_pressure!(::BoundarySPHSystem, ...)` is empty, + # so no pressure is updated in the previous update steps. update_pressure!(boundary_model, system, v, u, v_ode, u_ode, semi) return system diff --git a/src/schemes/fluid/weakly_compressible_sph/rhs.jl b/src/schemes/fluid/weakly_compressible_sph/rhs.jl index ecc6975d94..af088e0a1b 100644 --- a/src/schemes/fluid/weakly_compressible_sph/rhs.jl +++ b/src/schemes/fluid/weakly_compressible_sph/rhs.jl @@ -5,7 +5,7 @@ function interact!(dv, v_particle_system, u_particle_system, v_neighbor_system, u_neighbor_system, particle_system::WeaklyCompressibleSPHSystem, neighbor_system, semi) - (; density_calculator, state_equation, correction) = particle_system + (; density_calculator, correction) = particle_system sound_speed = system_sound_speed(particle_system) diff --git a/src/schemes/fluid/weakly_compressible_sph/system.jl b/src/schemes/fluid/weakly_compressible_sph/system.jl index e8de016119..bc9d541ea8 100644 --- a/src/schemes/fluid/weakly_compressible_sph/system.jl +++ b/src/schemes/fluid/weakly_compressible_sph/system.jl @@ -295,14 +295,16 @@ end function update_pressure!(system::WeaklyCompressibleSPHSystem, v, u, v_ode, u_ode, semi, t) (; density_calculator, correction, surface_normal_method, surface_tension) = system - compute_correction_values!(system, correction, u, v_ode, u_ode, semi) + compute_pressure!(system, v, semi) + # These are only computed when using corrections + compute_correction_values!(system, correction, u, v_ode, u_ode, semi) compute_gradient_correction_matrix!(correction, system, u, v_ode, u_ode, semi) - # `kernel_correct_density!` only performed for `SummationDensity` kernel_correct_density!(system, v, u, v_ode, u_ode, semi, correction, density_calculator) - compute_pressure!(system, v, semi) + + # 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_delta_function!(system, surface_tension, semi) return system From b3402ff563fbe33f00cf84209e75ba781612c7cf Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Fri, 1 Aug 2025 18:42:18 +0200 Subject: [PATCH 15/21] Fix boundary bug --- src/schemes/boundary/system.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/schemes/boundary/system.jl b/src/schemes/boundary/system.jl index 77607118f3..3ba368146d 100644 --- a/src/schemes/boundary/system.jl +++ b/src/schemes/boundary/system.jl @@ -359,7 +359,8 @@ end # This update depends on the computed quantities of the fluid system and therefore # has to be in `update_boundary_interpolation!` after `update_quantities!`. -function update_boundary_interpolation!(system::BoundarySPHSystem, v, u, v_ode, u_ode, semi) +function update_boundary_interpolation!(system::BoundarySPHSystem, v, u, v_ode, u_ode, + semi, t) (; boundary_model) = system # Note that `update_pressure!(::BoundarySPHSystem, ...)` is empty, From 5c6c8e5ff4d76c40f483a01aa129717b4a5b45d8 Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Sat, 2 Aug 2025 11:36:10 +0200 Subject: [PATCH 16/21] Make N-Body system type-stable --- examples/n_body/n_body_system.jl | 8 ++++---- src/general/semidiscretization.jl | 14 ++++++++++++-- 2 files changed, 16 insertions(+), 6 deletions(-) diff --git a/examples/n_body/n_body_system.jl b/examples/n_body/n_body_system.jl index 06e760fbe2..9493f2333f 100644 --- a/examples/n_body/n_body_system.jl +++ b/examples/n_body/n_body_system.jl @@ -1,8 +1,8 @@ using TrixiParticles using LinearAlgebra -struct NBodySystem{NDIMS, ELTYPE <: Real} <: TrixiParticles.System{NDIMS} - initial_condition :: InitialCondition{ELTYPE} +struct NBodySystem{NDIMS, ELTYPE <: Real, IC} <: TrixiParticles.System{NDIMS} + initial_condition :: IC mass :: Array{ELTYPE, 1} # [particle] G :: ELTYPE buffer :: Nothing @@ -11,13 +11,13 @@ struct NBodySystem{NDIMS, ELTYPE <: Real} <: TrixiParticles.System{NDIMS} mass = copy(initial_condition.mass) new{size(initial_condition.coordinates, 1), - eltype(mass)}(initial_condition, mass, G, nothing) + eltype(mass), typeof(initial_condition)}(initial_condition, mass, G, nothing) end end TrixiParticles.timer_name(::NBodySystem) = "nbody" -@inline Base.eltype(system::NBodySystem) = eltype(system.initial_condition.coordinates) +@inline Base.eltype(system::NBodySystem{NDIMS, ELTYPE}) where {NDIMS, ELTYPE} = ELTYPE function TrixiParticles.write_u0!(u0, system::NBodySystem) u0 .= system.initial_condition.coordinates diff --git a/src/general/semidiscretization.jl b/src/general/semidiscretization.jl index abe0d087f1..fa14c7d859 100644 --- a/src/general/semidiscretization.jl +++ b/src/general/semidiscretization.jl @@ -464,6 +464,18 @@ function drift!(du_ode, v_ode, u_ode, semi, t) end @inline function add_velocity!(du, v, particle, system) + # Generic fallback for all systems that don't define this function + for i in 1:ndims(system) + @inbounds du[i, particle] = v[i, particle] + end + + return du +end + +# Boundary systems don't integrate the particle positions +@inline add_velocity!(du, v, particle, system::BoundarySystem) = du + +@inline function add_velocity!(du, v, particle, system::FluidSystem) # This is zero unless a transport velocity is used delta_v_ = delta_v(system, particle) @@ -474,8 +486,6 @@ end return du end -@inline add_velocity!(du, v, particle, system::BoundarySPHSystem) = du - function kick!(dv_ode, v_ode, u_ode, semi, t) @trixi_timeit timer() "kick!" begin # Check that the `UpdateCallback` is used if required From ca0a2f6d3e131e8fdb389935cd46fe57cf40401e Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Sat, 2 Aug 2025 11:39:00 +0200 Subject: [PATCH 17/21] Fix packing --- src/callbacks/update.jl | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/src/callbacks/update.jl b/src/callbacks/update.jl index c3d5c59298..b989852203 100644 --- a/src/callbacks/update.jl +++ b/src/callbacks/update.jl @@ -89,6 +89,11 @@ function (update_callback!::UpdateCallback)(integrator) update_particle_packing(system, v_ode, u_ode, semi, integrator) end + # This is only used by the particle packing system and should be removed in the future + @trixi_timeit timer() "update TVF" foreach_system(semi) do system + update_transport_velocity!(system, v_ode, semi) + end + # Tell OrdinaryDiffEq that `u` has been modified u_modified!(integrator, true) From ea63b5fbd77964f69f0c5594b7bc49ba5317c072 Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Sat, 2 Aug 2025 14:40:36 +0200 Subject: [PATCH 18/21] Fix update callback --- src/preprocessing/particle_packing/system.jl | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/preprocessing/particle_packing/system.jl b/src/preprocessing/particle_packing/system.jl index 0e149178ea..fb97501e4c 100644 --- a/src/preprocessing/particle_packing/system.jl +++ b/src/preprocessing/particle_packing/system.jl @@ -390,6 +390,10 @@ end return system end +@inline function update_transport_velocity!(system, v_ode, semi) + return system +end + # Skip for fixed systems @inline add_velocity!(du, v, particle, system::ParticlePackingSystem{<:Any, true}) = du From 90865e00fb638f504f5b3e0a5952f568fb3d752e Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Sun, 3 Aug 2025 19:37:32 +0200 Subject: [PATCH 19/21] Update src/schemes/fluid/transport_velocity.jl Co-authored-by: Sven Berger --- src/schemes/fluid/transport_velocity.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/schemes/fluid/transport_velocity.jl b/src/schemes/fluid/transport_velocity.jl index bdfa0af8f1..bbf40bde2c 100644 --- a/src/schemes/fluid/transport_velocity.jl +++ b/src/schemes/fluid/transport_velocity.jl @@ -108,7 +108,7 @@ function update_tvf!(system, transport_velocity::TransportVelocityAdami, v, u, v semi; points=each_moving_particle(system)) do particle, neighbor, pos_diff, distance - m_a = @inbounds hydrodynamic_mass(neighbor_system, neighbor) + m_a = @inbounds hydrodynamic_mass(system, particle) m_b = @inbounds hydrodynamic_mass(neighbor_system, neighbor) rho_a = @inbounds current_density(v, system, particle) From d084d6caa02b5fc96661cfdc32b676a6b2f1113d Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Sun, 17 Aug 2025 15:58:52 +0300 Subject: [PATCH 20/21] Fix open boundaries --- src/general/semidiscretization.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/general/semidiscretization.jl b/src/general/semidiscretization.jl index fa14c7d859..cdc8cd34ea 100644 --- a/src/general/semidiscretization.jl +++ b/src/general/semidiscretization.jl @@ -472,8 +472,8 @@ end return du end -# Boundary systems don't integrate the particle positions -@inline add_velocity!(du, v, particle, system::BoundarySystem) = du +# Solid wall boundary system doesn't integrate the particle positions +@inline add_velocity!(du, v, particle, system::BoundarySPHSystem) = du @inline function add_velocity!(du, v, particle, system::FluidSystem) # This is zero unless a transport velocity is used From 92e9e7ac0e55f05073217efe9c5d8aa5cd194dcf Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Sun, 17 Aug 2025 16:03:44 +0300 Subject: [PATCH 21/21] Fix open boundaries --- src/schemes/boundary/open_boundary/system.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/schemes/boundary/open_boundary/system.jl b/src/schemes/boundary/open_boundary/system.jl index e820bd7aef..1d35bf726a 100644 --- a/src/schemes/boundary/open_boundary/system.jl +++ b/src/schemes/boundary/open_boundary/system.jl @@ -254,7 +254,7 @@ end end function update_boundary_interpolation!(system::OpenBoundarySPHSystem, v, u, v_ode, u_ode, - semi) + semi, t) update_boundary_model!(system, system.boundary_model, v, u, v_ode, u_ode, semi, t) end