From d4f676420d5317e4ede269d350f59ffb26f8465f Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Sat, 16 Aug 2025 19:57:23 +0300 Subject: [PATCH 01/13] Combine PST and TVF into a unified framework --- src/callbacks/callbacks.jl | 1 - src/callbacks/particle_shifting.jl | 148 -------------- src/callbacks/update.jl | 42 ++-- src/schemes/fluid/fluid.jl | 2 +- ...ort_velocity.jl => shifting_techniques.jl} | 186 ++++++++++++++---- .../fluid/weakly_compressible_sph/system.jl | 24 +-- 6 files changed, 186 insertions(+), 217 deletions(-) delete mode 100644 src/callbacks/particle_shifting.jl rename src/schemes/fluid/{transport_velocity.jl => shifting_techniques.jl} (52%) diff --git a/src/callbacks/callbacks.jl b/src/callbacks/callbacks.jl index 062a16e3a9..e9eb048c3f 100644 --- a/src/callbacks/callbacks.jl +++ b/src/callbacks/callbacks.jl @@ -32,4 +32,3 @@ include("post_process.jl") include("stepsize.jl") include("update.jl") include("steady_state_reached.jl") -include("particle_shifting.jl") diff --git a/src/callbacks/particle_shifting.jl b/src/callbacks/particle_shifting.jl deleted file mode 100644 index 921370f3b1..0000000000 --- a/src/callbacks/particle_shifting.jl +++ /dev/null @@ -1,148 +0,0 @@ -@doc raw""" - ParticleShiftingCallback() - -Callback to apply the Particle Shifting Technique by [Sun et al. (2017)](@cite Sun2017). -Following the original paper, the callback is applied in every time step and not -in every stage of a multi-stage time integration method to reduce the computational -cost and improve the stability of the scheme. - -See [Callbacks](@ref Callbacks) for more information on how to use this callback. -See [Particle Shifting Technique](@ref shifting) for more information on the method itself. - -!!! warning - The Particle Shifting Technique needs to be disabled close to the free surface - and therefore requires a free surface detection method. This is not yet implemented. - **This callback cannot be used in a free surface simulation.** -""" -function ParticleShiftingCallback() - # The first one is the `condition`, the second the `affect!` - return DiscreteCallback((particle_shifting_condition), particle_shifting!, - save_positions=(false, false)) -end - -# `condition` -function particle_shifting_condition(u, t, integrator) - return true -end - -# `affect!` -function particle_shifting!(integrator) - t = integrator.t - semi = integrator.p - v_ode, u_ode = integrator.u.x - dt = integrator.dt - # Internal cache vector, which is safe to use as temporary array - vu_cache = first(get_tmp_cache(integrator)) - - @trixi_timeit timer() "particle shifting callback" begin - # 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. - @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) - end - - @trixi_timeit timer() "particle shifting" foreach_system(semi) do system - u = wrap_u(u_ode, system, semi) - v = wrap_v(v_ode, system, semi) - particle_shifting!(u, v, system, v_ode, u_ode, semi, vu_cache, dt) - end - end - - # Tell OrdinaryDiffEq that `u` has been modified - u_modified!(integrator, true) - - return integrator -end - -function particle_shifting!(u, v, system, v_ode, u_ode, semi, u_cache, dt) - return u -end - -function particle_shifting!(u, v, system::FluidSystem, v_ode, u_ode, semi, - vu_cache, dt) - # Wrap the cache vector to an NDIMS x NPARTICLES matrix. - # We need this buffer because we cannot safely update `u` while iterating over it. - _, u_cache = vu_cache.x - delta_r = wrap_u(u_cache, system, semi) - set_zero!(delta_r) - - # This has similar performance to `maximum(..., eachparticle(system))`, - # but is GPU-compatible. - v_max = maximum(x -> sqrt(dot(x, x)), - reinterpret(reshape, SVector{ndims(system), eltype(v)}, - current_velocity(v, system))) - - # TODO this needs to be adapted to multi-resolution. - # Section 3.2 explains what else needs to be changed. - dx = particle_spacing(system, 1) - Wdx = smoothing_kernel(system, dx, 1) - h = smoothing_length(system, 1) - - foreach_system(semi) do neighbor_system - u_neighbor = wrap_u(u_ode, neighbor_system, semi) - v_neighbor = wrap_v(v_ode, neighbor_system, semi) - - system_coords = current_coordinates(u, system) - neighbor_coords = current_coordinates(u_neighbor, neighbor_system) - - foreach_point_neighbor(system, neighbor_system, system_coords, neighbor_coords, - semi; - points=each_moving_particle(system)) do particle, neighbor, - pos_diff, distance - m_b = hydrodynamic_mass(neighbor_system, neighbor) - rho_a = current_density(v, system, particle) - rho_b = current_density(v_neighbor, neighbor_system, neighbor) - - kernel = smoothing_kernel(system, distance, particle) - grad_kernel = smoothing_kernel_grad(system, pos_diff, distance, particle) - - # According to p. 29 below Eq. 9 - R = 2 // 10 - n = 4 - - # Eq. 7 in Sun et al. (2017). - # According to the paper, CFL * Ma can be rewritten as Δt * v_max / h - # (see p. 29, right above Eq. 9), but this does not work when scaling h. - # When setting CFL * Ma = Δt * v_max / (2 * Δx), PST works as expected - # for both small and large smoothing length factors. - # We need to scale - # - quadratically with the smoothing length, - # - linearly with the particle spacing, - # - linearly with the time step. - # See https://github.com/trixi-framework/TrixiParticles.jl/pull/834. - delta_r_ = -dt * v_max * (2 * h)^2 / (2 * dx) * (1 + R * (kernel / Wdx)^n) * - m_b / (rho_a + rho_b) * grad_kernel - - # Write into the buffer - for i in eachindex(delta_r_) - @inbounds delta_r[i, particle] += delta_r_[i] - end - end - end - - # Add δ_r from the buffer to the current coordinates - @threaded semi for particle in eachparticle(system) - for i in axes(delta_r, 1) - @inbounds u[i, particle] += delta_r[i, particle] - end - end - - return u -end - -function Base.show(io::IO, cb::DiscreteCallback{<:Any, typeof(particle_shifting!)}) - @nospecialize cb # reduce precompilation time - print(io, "ParticleShiftingCallback()") -end - -function Base.show(io::IO, ::MIME"text/plain", - cb::DiscreteCallback{<:Any, typeof(particle_shifting!)}) - @nospecialize cb # reduce precompilation time - - if get(io, :compact, false) - show(io, cb) - else - summary_box(io, "ParticleShiftingCallback") - end -end diff --git a/src/callbacks/update.jl b/src/callbacks/update.jl index 77496a46c5..1d348f76c1 100644 --- a/src/callbacks/update.jl +++ b/src/callbacks/update.jl @@ -74,22 +74,32 @@ function (update_callback!::UpdateCallback)(integrator) semi = integrator.p v_ode, u_ode = integrator.u.x - # 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 open boundaries first, since particles might be activated or deactivated - @trixi_timeit timer() "update open boundary" foreach_system(semi) do system - update_open_boundary_eachstep!(system, v_ode, u_ode, semi, t) - end - - @trixi_timeit timer() "update particle packing" foreach_system(semi) do system - 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) + @trixi_timeit timer() "update callback" begin + # 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. + @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) + end + + # Update open boundaries first, since particles might be activated or deactivated + @trixi_timeit timer() "update open boundary" foreach_system(semi) do system + update_open_boundary_eachstep!(system, v_ode, u_ode, semi, t) + end + + @trixi_timeit timer() "update particle packing" foreach_system(semi) do system + 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 + + @trixi_timeit timer() "particle shifting" foreach_system(semi) do system + particle_shifting!(u_ode, shifting_technique(system), system, semi, + integrator.dt) + end end # Tell OrdinaryDiffEq that `u` has been modified diff --git a/src/schemes/fluid/fluid.jl b/src/schemes/fluid/fluid.jl index d3e6adf971..59703f9b6d 100644 --- a/src/schemes/fluid/fluid.jl +++ b/src/schemes/fluid/fluid.jl @@ -193,7 +193,7 @@ end include("pressure_acceleration.jl") include("viscosity.jl") -include("transport_velocity.jl") +include("shifting_techniques.jl") include("surface_tension.jl") include("surface_normal_sph.jl") include("weakly_compressible_sph/weakly_compressible_sph.jl") diff --git a/src/schemes/fluid/transport_velocity.jl b/src/schemes/fluid/shifting_techniques.jl similarity index 52% rename from src/schemes/fluid/transport_velocity.jl rename to src/schemes/fluid/shifting_techniques.jl index fce7375dd5..4ebbe01758 100644 --- a/src/schemes/fluid/transport_velocity.jl +++ b/src/schemes/fluid/shifting_techniques.jl @@ -1,56 +1,168 @@ -""" - TransportVelocityAdami(background_pressure::Real) - -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 -- `background_pressure`: Background pressure. Suggested is a background pressure which is - on the order of the reference pressure. -""" -struct TransportVelocityAdami{T <: Real} - background_pressure::T -end +abstract type AbstractShiftingTechnique end -# No TVF for a system by default -@inline transport_velocity(system) = nothing +# No shifting for a system by default +@inline shifting_technique(system) = nothing -create_cache_tvf(initial_condition, ::Nothing) = (;) +create_cache_shifting(initial_condition, ::Nothing) = (;) -function create_cache_tvf(initial_condition, ::TransportVelocityAdami) +function create_cache_shifting(initial_condition, ::AbstractShiftingTechnique) delta_v = zeros(eltype(initial_condition), ndims(initial_condition), nparticles(initial_condition)) return (; delta_v) end -# `δv` is the correction to the particle velocity due to the TVF. +# `δv` is the correction to the particle velocity due to the shifting. # Particles are advected with the velocity `v + δv`. @propagate_inbounds function delta_v(system, particle) - return delta_v(system, transport_velocity(system), particle) + return delta_v(system, shifting_technique(system), particle) +end + +# Zero when no shifting is used +@inline function delta_v(system, shifting, particle) + return zero(SVector{ndims(system), eltype(system)}) end -@propagate_inbounds function delta_v(system, ::TransportVelocityAdami, particle) +@propagate_inbounds function delta_v(system, ::AbstractShiftingTechnique, particle) return extract_svector(system.cache.delta_v, system, particle) end -# Zero when no TVF is used -@inline function delta_v(system, transport_velocity, particle) - return zero(SVector{ndims(system), eltype(system)}) +function update_shifting!(system, shifting, v, u, v_ode, u_ode, semi, t) + return system end -@inline function dv_transport_velocity(::Nothing, system, neighbor_system, - particle, neighbor, v_system, v_neighbor_system, - m_a, m_b, rho_a, rho_b, pos_diff, distance, - grad_kernel, correction) +# Additional term in the momentum equation due to the shifting technique +@inline function dv_shifting(shifting, system, neighbor_system, + 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, - particle, neighbor, v_system, v_neighbor_system, - m_a, m_b, rho_a, rho_b, pos_diff, distance, - grad_kernel, correction) +@doc raw""" + ParticleShiftingTechnique() + +Particle Shifting Technique by [Sun et al. (2017)](@cite Sun2017). +Following the original paper, the callback is applied in every time step and not +in every stage of a multi-stage time integration method to reduce the computational +cost and improve the stability of the scheme. + +See [Particle Shifting Technique](@ref shifting) for more information on the method. + +!!! warning + The Particle Shifting Technique needs to be disabled close to the free surface + and therefore requires a free surface detection method. This is not yet implemented. + **This technique cannot be used in a free surface simulation.** +""" +struct ParticleShiftingTechnique end + +# Zero because PST is applied in a callback +@inline function delta_v(system, ::ParticleShiftingTechnique, particle) + return zero(SVector{ndims(system), eltype(system)}) +end + +# This is called from the update callback +particle_shifting!(u_ode, shifting, system, semi, dt) = u_ode + +function particle_shifting!(u_ode, ::ParticleShiftingTechnique, system, semi, dt) + (; cache) = system + (; delta_v) = cache + + u = wrap_u(u_ode, system, semi) + + # Add δr from the cache to the current coordinates + @threaded semi for particle in eachparticle(system) + for i in axes(delta_v, 1) + @inbounds u[i, particle] += dt * delta_v[i, particle] + end + end + + return u +end + +function update_shifting!(system, ::ParticleShiftingTechnique, + v, u, v_ode, u_ode, semi, t) + (; cache) = system + (; delta_v) = cache + + set_zero!(delta_v) + + # This has similar performance to `maximum(..., eachparticle(system))`, + # but is GPU-compatible. + v_max = maximum(x -> sqrt(dot(x, x)), + reinterpret(reshape, SVector{ndims(system), eltype(v)}, + current_velocity(v, system))) + + # TODO this needs to be adapted to multi-resolution. + # Section 3.2 explains what else needs to be changed. + dx = particle_spacing(system, 1) + Wdx = smoothing_kernel(system, dx, 1) + h = smoothing_length(system, 1) + + foreach_system(semi) do neighbor_system + u_neighbor = wrap_u(u_ode, neighbor_system, semi) + v_neighbor = wrap_v(v_ode, neighbor_system, semi) + + system_coords = current_coordinates(u, system) + neighbor_coords = current_coordinates(u_neighbor, neighbor_system) + + foreach_point_neighbor(system, neighbor_system, system_coords, neighbor_coords, + semi; + points=each_moving_particle(system)) do particle, neighbor, + pos_diff, distance + m_b = hydrodynamic_mass(neighbor_system, neighbor) + rho_a = current_density(v, system, particle) + rho_b = current_density(v_neighbor, neighbor_system, neighbor) + + kernel = smoothing_kernel(system, distance, particle) + grad_kernel = smoothing_kernel_grad(system, pos_diff, distance, particle) + + # According to p. 29 below Eq. 9 + R = 2 // 10 + n = 4 + + # Eq. 7 in Sun et al. (2017). + # According to the paper, CFL * Ma can be rewritten as Δt * v_max / h + # (see p. 29, right above Eq. 9), but this does not work when scaling h. + # When setting CFL * Ma = Δt * v_max / (2 * Δx), PST works as expected + # for both small and large smoothing length factors. + # We need to scale + # - quadratically with the smoothing length, + # - linearly with the particle spacing, + # - linearly with the time step. + # See https://github.com/trixi-framework/TrixiParticles.jl/pull/834. + delta_v_ = -v_max * (2 * h)^2 / (2 * dx) * (1 + R * (kernel / Wdx)^n) * + m_b / (rho_a + rho_b) * grad_kernel + + # Write into the buffer + for i in eachindex(delta_v_) + @inbounds delta_v[i, particle] += delta_v_[i] + end + end + end + + return system +end + +""" + TransportVelocityAdami(background_pressure::Real) + +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 +- `background_pressure`: Background pressure. Suggested is a background pressure which is + on the order of the reference pressure. +""" +struct TransportVelocityAdami{T <: Real} + background_pressure::T +end + +@inline function dv_shifting(::TransportVelocityAdami, system, neighbor_system, + 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) @@ -72,15 +184,11 @@ end distance, grad_kernel, correction) end -function update_tvf!(system, transport_velocity, v, u, v_ode, u_ode, semi, t) - return system -end - -function update_tvf!(system, transport_velocity::TransportVelocityAdami, v, u, v_ode, - u_ode, semi, t) +function update_shifting!(system, shifting::TransportVelocityAdami, v, u, v_ode, + u_ode, semi, t) (; cache, correction) = system (; delta_v) = cache - (; background_pressure) = transport_velocity + (; background_pressure) = shifting sound_speed = system_sound_speed(system) diff --git a/src/schemes/fluid/weakly_compressible_sph/system.jl b/src/schemes/fluid/weakly_compressible_sph/system.jl index a5f56499b8..884dd24105 100644 --- a/src/schemes/fluid/weakly_compressible_sph/system.jl +++ b/src/schemes/fluid/weakly_compressible_sph/system.jl @@ -5,7 +5,7 @@ acceleration=ntuple(_ -> 0.0, NDIMS), viscosity=nothing, density_diffusion=nothing, pressure_acceleration=nothing, - transport_velocity=nothing, + shifting_technique=nothing, buffer_size=nothing, correction=nothing, source_terms=nothing, surface_tension=nothing, surface_normal_method=nothing, @@ -36,8 +36,8 @@ See [Weakly Compressible SPH](@ref wcsph) for more details on the method. density calculator and the correction method. To use [Tensile Instability Control](@ref tic), pass [`tensile_instability_control`](@ref) here. -- `transport_velocity`: [Transport Velocity Formulation (TVF)](@ref transport_velocity_formulation). - Default is no TVF. +- `shifting_technique`: [Shifting Correction](@ref shifting_technique). + Default is no shifting. - `buffer_size`: Number of buffer particles. This is needed when simulating with [`OpenBoundarySPHSystem`](@ref). - `correction`: Correction method used for this system. (default: no correction, see [Corrections](@ref corrections)) @@ -59,7 +59,7 @@ See [Weakly Compressible SPH](@ref wcsph) for more details on the method. - `color_value`: The value used to initialize the color of particles in the system. """ struct WeaklyCompressibleSPHSystem{NDIMS, ELTYPE <: Real, IC, MA, P, DC, SE, K, V, DD, COR, - PF, TV, ST, B, SRFT, SRFN, PR, C} <: FluidSystem{NDIMS} + PF, SC, ST, B, SRFT, SRFN, PR, C} <: FluidSystem{NDIMS} initial_condition :: IC mass :: MA # Array{ELTYPE, 1} pressure :: P # Array{ELTYPE, 1} @@ -71,7 +71,7 @@ struct WeaklyCompressibleSPHSystem{NDIMS, ELTYPE <: Real, IC, MA, P, DC, SE, K, density_diffusion :: DD correction :: COR pressure_acceleration_formulation :: PF - transport_velocity :: TV + shifting_technique :: SC source_terms :: ST surface_tension :: SRFT surface_normal_method :: SRFN @@ -89,7 +89,7 @@ function WeaklyCompressibleSPHSystem(initial_condition, ndims(smoothing_kernel)), viscosity=nothing, density_diffusion=nothing, pressure_acceleration=nothing, - transport_velocity=nothing, + shifting_technique=nothing, buffer_size=nothing, correction=nothing, source_terms=nothing, surface_tension=nothing, surface_normal_method=nothing, @@ -147,7 +147,7 @@ function WeaklyCompressibleSPHSystem(initial_condition, n_particles)..., create_cache_refinement(initial_condition, particle_refinement, smoothing_length)..., - create_cache_tvf(initial_condition, transport_velocity)..., + create_cache_shifting(initial_condition, shifting_technique)..., color=Int(color_value)) # If the `reference_density_spacing` is set calculate the `ideal_neighbor_count` @@ -162,7 +162,7 @@ function WeaklyCompressibleSPHSystem(initial_condition, density_calculator, state_equation, smoothing_kernel, acceleration_, viscosity, density_diffusion, correction, pressure_acceleration, - transport_velocity, source_terms, surface_tension, + shifting_technique, source_terms, surface_tension, surface_normal_method, buffer, particle_refinement, cache) end @@ -177,6 +177,7 @@ function Base.show(io::IO, system::WeaklyCompressibleSPHSystem) print(io, ", ", system.smoothing_kernel) print(io, ", ", system.viscosity) print(io, ", ", system.density_diffusion) + print(io, ", ", system.shifting_technique) print(io, ", ", system.surface_tension) print(io, ", ", system.surface_normal_method) if system.surface_normal_method isa ColorfieldSurfaceNormal @@ -207,9 +208,8 @@ function Base.show(io::IO, ::MIME"text/plain", system::WeaklyCompressibleSPHSyst summary_line(io, "state equation", system.state_equation |> typeof |> nameof) summary_line(io, "smoothing kernel", system.smoothing_kernel |> typeof |> nameof) summary_line(io, "viscosity", system.viscosity) - summary_line(io, "tansport velocity formulation", - system.transport_velocity |> typeof |> nameof) summary_line(io, "density diffusion", system.density_diffusion) + summary_line(io, "shifting correction", system.shifting_technique) summary_line(io, "surface tension", system.surface_tension) summary_line(io, "surface normal method", system.surface_normal_method) if system.surface_normal_method isa ColorfieldSurfaceNormal @@ -278,7 +278,7 @@ end @inline system_sound_speed(system::WeaklyCompressibleSPHSystem) = system.state_equation.sound_speed -@inline transport_velocity(system::WeaklyCompressibleSPHSystem) = system.transport_velocity +@inline shifting_technique(system::WeaklyCompressibleSPHSystem) = system.shifting_technique function update_quantities!(system::WeaklyCompressibleSPHSystem, v, u, v_ode, u_ode, semi, t) @@ -316,7 +316,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) - update_tvf!(system, transport_velocity(system), v, u, v_ode, u_ode, semi, t) + update_shifting!(system, shifting_technique(system), v, u, v_ode, u_ode, semi, t) end function kernel_correct_density!(system::WeaklyCompressibleSPHSystem, v, u, v_ode, u_ode, From 62805389a35f81e3b06fc62f8c2fc257625f2a4c Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Sat, 16 Aug 2025 20:21:29 +0300 Subject: [PATCH 02/13] Require update callback for PST --- src/general/system.jl | 3 --- src/schemes/fluid/shifting_techniques.jl | 6 ++++++ 2 files changed, 6 insertions(+), 3 deletions(-) diff --git a/src/general/system.jl b/src/general/system.jl index 929fb8b425..5b4d5aecb8 100644 --- a/src/general/system.jl +++ b/src/general/system.jl @@ -147,9 +147,6 @@ function update_final!(system, v, u, v_ode, u_ode, semi, t) return system end -# Only for systems requiring the use of the `UpdateCallback` -@inline requires_update_callback(system) = false - @inline initial_smoothing_length(system) = smoothing_length(system, nothing) @inline function smoothing_length(system, particle) diff --git a/src/schemes/fluid/shifting_techniques.jl b/src/schemes/fluid/shifting_techniques.jl index 4ebbe01758..9e3892068a 100644 --- a/src/schemes/fluid/shifting_techniques.jl +++ b/src/schemes/fluid/shifting_techniques.jl @@ -3,6 +3,12 @@ abstract type AbstractShiftingTechnique end # No shifting for a system by default @inline shifting_technique(system) = nothing +# WARNING: Be careful if defining this function for a specific system type. +# The version for a specific system type will override this generic version. +requires_update_callback(system) = requires_update_callback(shifting_technique(system)) +requires_update_callback(::Nothing) = false +requires_update_callback(::AbstractShiftingTechnique) = true + create_cache_shifting(initial_condition, ::Nothing) = (;) function create_cache_shifting(initial_condition, ::AbstractShiftingTechnique) From 1f882f6ac7251a3e79203ac13032fc239dd292f4 Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Sun, 17 Aug 2025 14:08:06 +0300 Subject: [PATCH 03/13] Fix WCSPH --- src/schemes/fluid/shifting_techniques.jl | 4 ++-- src/schemes/fluid/weakly_compressible_sph/rhs.jl | 11 +++++------ 2 files changed, 7 insertions(+), 8 deletions(-) diff --git a/src/schemes/fluid/shifting_techniques.jl b/src/schemes/fluid/shifting_techniques.jl index 9e3892068a..e2af998397 100644 --- a/src/schemes/fluid/shifting_techniques.jl +++ b/src/schemes/fluid/shifting_techniques.jl @@ -60,7 +60,7 @@ See [Particle Shifting Technique](@ref shifting) for more information on the met and therefore requires a free surface detection method. This is not yet implemented. **This technique cannot be used in a free surface simulation.** """ -struct ParticleShiftingTechnique end +struct ParticleShiftingTechnique <: AbstractShiftingTechnique end # Zero because PST is applied in a callback @inline function delta_v(system, ::ParticleShiftingTechnique, particle) @@ -161,7 +161,7 @@ See [TVF](@ref transport_velocity_formulation) for more details of the method. - `background_pressure`: Background pressure. Suggested is a background pressure which is on the order of the reference pressure. """ -struct TransportVelocityAdami{T <: Real} +struct TransportVelocityAdami{T <: Real} <: AbstractShiftingTechnique background_pressure::T end diff --git a/src/schemes/fluid/weakly_compressible_sph/rhs.jl b/src/schemes/fluid/weakly_compressible_sph/rhs.jl index af088e0a1b..360c5ec22d 100644 --- a/src/schemes/fluid/weakly_compressible_sph/rhs.jl +++ b/src/schemes/fluid/weakly_compressible_sph/rhs.jl @@ -71,12 +71,11 @@ function interact!(dv, v_particle_system, u_particle_system, grad_kernel) # Add convection term (only when using `TransportVelocityAdami`) - 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_tvf = dv_shifting(shifting_technique(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, From 8b150bb0c6a86572cbcb97f7ab86652998f53697 Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Sun, 17 Aug 2025 15:05:26 +0300 Subject: [PATCH 04/13] Update PST only in callback --- src/callbacks/update.jl | 4 +- src/schemes/fluid/shifting_techniques.jl | 50 ++++++++++++------- .../fluid/weakly_compressible_sph/system.jl | 4 +- 3 files changed, 36 insertions(+), 22 deletions(-) diff --git a/src/callbacks/update.jl b/src/callbacks/update.jl index 1d348f76c1..8eecb8380a 100644 --- a/src/callbacks/update.jl +++ b/src/callbacks/update.jl @@ -97,8 +97,8 @@ function (update_callback!::UpdateCallback)(integrator) end @trixi_timeit timer() "particle shifting" foreach_system(semi) do system - particle_shifting!(u_ode, shifting_technique(system), system, semi, - integrator.dt) + particle_shifting_from_callback!(u_ode, shifting_technique(system), system, + v_ode, semi, integrator.dt) end end diff --git a/src/schemes/fluid/shifting_techniques.jl b/src/schemes/fluid/shifting_techniques.jl index e2af998397..13c78a4f9c 100644 --- a/src/schemes/fluid/shifting_techniques.jl +++ b/src/schemes/fluid/shifting_techniques.jl @@ -9,6 +9,9 @@ requires_update_callback(system) = requires_update_callback(shifting_technique(s requires_update_callback(::Nothing) = false requires_update_callback(::AbstractShiftingTechnique) = true +# This is called from the `UpdateCallback` +particle_shifting_from_callback!(u_ode, shifting, system, v_ode, semi, dt) = u_ode + create_cache_shifting(initial_condition, ::Nothing) = (;) function create_cache_shifting(initial_condition, ::AbstractShiftingTechnique) @@ -33,7 +36,7 @@ end return extract_svector(system.cache.delta_v, system, particle) end -function update_shifting!(system, shifting, v, u, v_ode, u_ode, semi, t) +function update_shifting!(system, shifting, v, u, v_ode, u_ode, semi) return system end @@ -67,27 +70,22 @@ struct ParticleShiftingTechnique <: AbstractShiftingTechnique end return zero(SVector{ndims(system), eltype(system)}) end -# This is called from the update callback -particle_shifting!(u_ode, shifting, system, semi, dt) = u_ode - -function particle_shifting!(u_ode, ::ParticleShiftingTechnique, system, semi, dt) - (; cache) = system - (; delta_v) = cache +function particle_shifting_from_callback!(u_ode, shifting::ParticleShiftingTechnique, + system, v_ode, semi, dt) + @trixi_timeit timer() "particle shifting" begin + v = wrap_v(v_ode, system, semi) + u = wrap_u(u_ode, system, semi) - u = wrap_u(u_ode, system, semi) + # Update the shifting velocity + update_shifting_from_callback!(system, shifting, v, u, v_ode, u_ode, semi) - # Add δr from the cache to the current coordinates - @threaded semi for particle in eachparticle(system) - for i in axes(delta_v, 1) - @inbounds u[i, particle] += dt * delta_v[i, particle] - end + # Update the particle positions with the shifting velocity + particle_shifting!(u_ode, shifting, system, semi, dt) end - - return u end -function update_shifting!(system, ::ParticleShiftingTechnique, - v, u, v_ode, u_ode, semi, t) +function update_shifting_from_callback!(system, ::ParticleShiftingTechnique, + v, u, v_ode, u_ode, semi) (; cache) = system (; delta_v) = cache @@ -150,6 +148,22 @@ function update_shifting!(system, ::ParticleShiftingTechnique, return system end +function particle_shifting!(u_ode, ::ParticleShiftingTechnique, system, semi, dt) + (; cache) = system + (; delta_v) = cache + + u = wrap_u(u_ode, system, semi) + + # Add δr from the cache to the current coordinates + @threaded semi for particle in eachparticle(system) + for i in axes(delta_v, 1) + @inbounds u[i, particle] += dt * delta_v[i, particle] + end + end + + return u +end + """ TransportVelocityAdami(background_pressure::Real) @@ -191,7 +205,7 @@ end end function update_shifting!(system, shifting::TransportVelocityAdami, v, u, v_ode, - u_ode, semi, t) + u_ode, semi) (; cache, correction) = system (; delta_v) = cache (; background_pressure) = shifting diff --git a/src/schemes/fluid/weakly_compressible_sph/system.jl b/src/schemes/fluid/weakly_compressible_sph/system.jl index 884dd24105..a6747822bc 100644 --- a/src/schemes/fluid/weakly_compressible_sph/system.jl +++ b/src/schemes/fluid/weakly_compressible_sph/system.jl @@ -71,7 +71,7 @@ struct WeaklyCompressibleSPHSystem{NDIMS, ELTYPE <: Real, IC, MA, P, DC, SE, K, density_diffusion :: DD correction :: COR pressure_acceleration_formulation :: PF - shifting_technique :: SC + shifting_technique :: SC source_terms :: ST surface_tension :: SRFT surface_normal_method :: SRFN @@ -316,7 +316,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) - update_shifting!(system, shifting_technique(system), v, u, v_ode, u_ode, semi, t) + update_shifting!(system, shifting_technique(system), v, u, v_ode, u_ode, semi) end function kernel_correct_density!(system::WeaklyCompressibleSPHSystem, v, u, v_ode, u_ode, From fa90ae975378b23cf2c18abf619b66a81265a945 Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Sun, 17 Aug 2025 15:22:59 +0300 Subject: [PATCH 05/13] Fix EDAC --- src/io/write_vtk.jl | 3 -- .../fluid/entropically_damped_sph/rhs.jl | 13 ++++----- .../fluid/entropically_damped_sph/system.jl | 28 +++++++++---------- .../fluid/weakly_compressible_sph/rhs.jl | 2 +- .../fluid/weakly_compressible_sph/system.jl | 7 +++-- test/systems/edac_system.jl | 4 +-- 6 files changed, 27 insertions(+), 30 deletions(-) diff --git a/src/io/write_vtk.jl b/src/io/write_vtk.jl index b93e25f324..1aca732098 100644 --- a/src/io/write_vtk.jl +++ b/src/io/write_vtk.jl @@ -349,9 +349,6 @@ function write2vtk!(vtk, v, u, t, system::FluidSystem; write_meta_data=true) else vtk["solver"] = "EDAC" vtk["sound_speed"] = system.sound_speed - vtk["background_pressure_TVF"] = system.transport_velocity isa Nothing ? - "-" : - system.transport_velocity.background_pressure end end diff --git a/src/schemes/fluid/entropically_damped_sph/rhs.jl b/src/schemes/fluid/entropically_damped_sph/rhs.jl index 08139e2ca6..a4d7d280a4 100644 --- a/src/schemes/fluid/entropically_damped_sph/rhs.jl +++ b/src/schemes/fluid/entropically_damped_sph/rhs.jl @@ -51,13 +51,12 @@ function interact!(dv, v_particle_system, u_particle_system, particle, neighbor, pos_diff, distance, sound_speed, m_a, m_b, rho_a, rho_b, grad_kernel) - # Add convection term (only when using `TransportVelocityAdami`) - 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) + # Extra terms in the momentum equation when using a shifting technique + dv_tvf = dv_shifting(shifting_technique(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, diff --git a/src/schemes/fluid/entropically_damped_sph/system.jl b/src/schemes/fluid/entropically_damped_sph/system.jl index 19e50c441a..89d8a923ff 100644 --- a/src/schemes/fluid/entropically_damped_sph/system.jl +++ b/src/schemes/fluid/entropically_damped_sph/system.jl @@ -3,7 +3,7 @@ smoothing_length, sound_speed; pressure_acceleration=inter_particle_averaged_pressure, density_calculator=SummationDensity(), - transport_velocity=nothing, + shifting_technique=nothing, alpha=0.5, viscosity=nothing, acceleration=ntuple(_ -> 0.0, NDIMS), surface_tension=nothing, surface_normal_method=nothing, buffer_size=nothing, @@ -30,10 +30,11 @@ See [Entropically Damped Artificial Compressibility for SPH](@ref edac) for more When set to `nothing`, the pressure acceleration formulation for the corresponding [density calculator](@ref density_calculator) is chosen. - `density_calculator`: [Density calculator](@ref density_calculator) (default: [`SummationDensity`](@ref)) -- `transport_velocity`: [Transport Velocity Formulation (TVF)](@ref transport_velocity_formulation). - Default is no TVF. +- `shifting_technique`: [Shifting technique](@ref shifting) or [transport velocity + formulation](@ref transport_velocity_formulation) to use + with this system. Default is no shifting. - `average_pressure_reduction`: Whether to subtract the average pressure of neighboring particles - from the local pressure (default: `true` when using TVF, `false` otherwise). + from the local pressure (default: `true` when using shifting, `false` otherwise). - `buffer_size`: Number of buffer particles. This is needed when simulating with [`OpenBoundarySPHSystem`](@ref). - `correction`: Correction method used for this system. (default: no correction, see [Corrections](@ref corrections)) @@ -67,7 +68,7 @@ struct EntropicallyDampedSPHSystem{NDIMS, ELTYPE <: Real, IC, M, DC, K, V, COR, acceleration :: SVector{NDIMS, ELTYPE} correction :: COR pressure_acceleration_formulation :: PF - transport_velocity :: TV + shifting_technique :: TV average_pressure_reduction :: AVGP source_terms :: ST surface_tension :: SRFT @@ -83,8 +84,8 @@ function EntropicallyDampedSPHSystem(initial_condition, smoothing_kernel, smoothing_length, sound_speed; pressure_acceleration=inter_particle_averaged_pressure, density_calculator=SummationDensity(), - transport_velocity=nothing, - average_pressure_reduction=(!isnothing(transport_velocity)), + shifting_technique=nothing, + average_pressure_reduction=(!isnothing(shifting_technique)), alpha=0.5, viscosity=nothing, acceleration=ntuple(_ -> 0.0, ndims(smoothing_kernel)), @@ -137,7 +138,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(initial_condition, transport_velocity)..., + create_cache_shifting(initial_condition, shifting_technique)..., create_cache_avg_pressure_reduction(initial_condition, avg_pressure_reduction)..., create_cache_surface_normal(surface_normal_method, ELTYPE, NDIMS, @@ -161,14 +162,14 @@ function EntropicallyDampedSPHSystem(initial_condition, smoothing_kernel, EntropicallyDampedSPHSystem{NDIMS, ELTYPE, typeof(initial_condition), typeof(mass), typeof(density_calculator), typeof(smoothing_kernel), typeof(viscosity), typeof(correction), - typeof(pressure_acceleration), typeof(transport_velocity), + typeof(pressure_acceleration), typeof(shifting_technique), typeof(avg_pressure_reduction), typeof(source_terms), typeof(surface_tension), typeof(surface_normal_method), typeof(buffer), Nothing, typeof(cache)}(initial_condition, mass, density_calculator, smoothing_kernel, sound_speed, viscosity, nu_edac, acceleration_, correction, - pressure_acceleration, transport_velocity, + pressure_acceleration, shifting_technique, avg_pressure_reduction, source_terms, surface_tension, surface_normal_method, buffer, @@ -218,8 +219,7 @@ function Base.show(io::IO, ::MIME"text/plain", system::EntropicallyDampedSPHSyst summary_line(io, "viscosity", system.viscosity |> typeof |> nameof) summary_line(io, "ν₍EDAC₎", "≈ $(round(system.nu_edac; digits=3))") summary_line(io, "smoothing kernel", system.smoothing_kernel |> typeof |> nameof) - summary_line(io, "tansport velocity formulation", - system.transport_velocity |> typeof |> nameof) + summary_line(io, "shifting technique", system.shifting_technique) summary_line(io, "average pressure reduction", typeof(system.average_pressure_reduction).parameters[1] ? "yes" : "no") summary_line(io, "acceleration", system.acceleration) @@ -271,7 +271,7 @@ end @inline system_sound_speed(system::EntropicallyDampedSPHSystem) = system.sound_speed -@inline transport_velocity(system::EntropicallyDampedSPHSystem) = system.transport_velocity +@inline shifting_technique(system::EntropicallyDampedSPHSystem) = system.shifting_technique @inline function average_pressure(system::EntropicallyDampedSPHSystem, particle) average_pressure(system, system.average_pressure_reduction, particle) @@ -321,7 +321,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.average_pressure_reduction, v_ode, u_ode, semi) - update_tvf!(system, transport_velocity(system), v, u, v_ode, u_ode, semi, t) + update_shifting!(system, shifting_technique(system), v, u, v_ode, u_ode, semi) end # No average pressure reduction is used diff --git a/src/schemes/fluid/weakly_compressible_sph/rhs.jl b/src/schemes/fluid/weakly_compressible_sph/rhs.jl index 360c5ec22d..b9cd27f079 100644 --- a/src/schemes/fluid/weakly_compressible_sph/rhs.jl +++ b/src/schemes/fluid/weakly_compressible_sph/rhs.jl @@ -70,7 +70,7 @@ 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`) + # Extra terms in the momentum equation when using a shifting technique dv_tvf = dv_shifting(shifting_technique(particle_system), particle_system, neighbor_system, particle, neighbor, v_particle_system, v_neighbor_system, diff --git a/src/schemes/fluid/weakly_compressible_sph/system.jl b/src/schemes/fluid/weakly_compressible_sph/system.jl index a6747822bc..dda2ead5f1 100644 --- a/src/schemes/fluid/weakly_compressible_sph/system.jl +++ b/src/schemes/fluid/weakly_compressible_sph/system.jl @@ -36,8 +36,9 @@ See [Weakly Compressible SPH](@ref wcsph) for more details on the method. density calculator and the correction method. To use [Tensile Instability Control](@ref tic), pass [`tensile_instability_control`](@ref) here. -- `shifting_technique`: [Shifting Correction](@ref shifting_technique). - Default is no shifting. +- `shifting_technique`: [Shifting technique](@ref shifting) or [transport velocity + formulation](@ref transport_velocity_formulation) to use + with this system. Default is no shifting. - `buffer_size`: Number of buffer particles. This is needed when simulating with [`OpenBoundarySPHSystem`](@ref). - `correction`: Correction method used for this system. (default: no correction, see [Corrections](@ref corrections)) @@ -209,7 +210,7 @@ function Base.show(io::IO, ::MIME"text/plain", system::WeaklyCompressibleSPHSyst summary_line(io, "smoothing kernel", system.smoothing_kernel |> typeof |> nameof) summary_line(io, "viscosity", system.viscosity) summary_line(io, "density diffusion", system.density_diffusion) - summary_line(io, "shifting correction", system.shifting_technique) + summary_line(io, "shifting technique", system.shifting_technique) summary_line(io, "surface tension", system.surface_tension) summary_line(io, "surface normal method", system.surface_normal_method) if system.surface_normal_method isa ColorfieldSurfaceNormal diff --git a/test/systems/edac_system.jl b/test/systems/edac_system.jl index 3d8cf9babb..cdebb91768 100644 --- a/test/systems/edac_system.jl +++ b/test/systems/edac_system.jl @@ -32,7 +32,7 @@ @test system.mass == mass @test system.smoothing_kernel == smoothing_kernel @test TrixiParticles.initial_smoothing_length(system) == smoothing_length - @test system.transport_velocity isa Nothing + @test system.shifting_technique isa Nothing @test system.viscosity === nothing @test system.nu_edac == (0.5 * smoothing_length * sound_speed) / 8 @test system.acceleration == [0.0 for _ in 1:NDIMS] @@ -87,7 +87,7 @@ @test system.mass == setup.mass @test system.smoothing_kernel == smoothing_kernel @test TrixiParticles.initial_smoothing_length(system) == smoothing_length - @test system.transport_velocity isa Nothing + @test system.shifting_technique isa Nothing @test system.viscosity === nothing @test system.nu_edac == (0.5 * smoothing_length * sound_speed) / 8 @test system.acceleration == [0.0 for _ in 1:NDIMS] From 9bdd2d6f276e8b083b208c8464f1c63363b462c2 Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Sun, 17 Aug 2025 16:25:33 +0300 Subject: [PATCH 06/13] Update docs --- docs/src/systems/entropically_damped_sph.md | 61 ------------------ docs/src/systems/weakly_compressible_sph.md | 70 ++++++++++++++++++++- 2 files changed, 68 insertions(+), 63 deletions(-) diff --git a/docs/src/systems/entropically_damped_sph.md b/docs/src/systems/entropically_damped_sph.md index 0db8d275e2..d035fb6fc1 100644 --- a/docs/src/systems/entropically_damped_sph.md +++ b/docs/src/systems/entropically_damped_sph.md @@ -45,64 +45,3 @@ is a good choice for a wide range of Reynolds numbers (0.0125 to 10000). Modules = [TrixiParticles] Pages = [joinpath("schemes", "fluid", "entropically_damped_sph", "system.jl")] ``` - -## [Transport Velocity Formulation (TVF)](@id transport_velocity_formulation) -Standard SPH suffers from problems like tensile instability or the creation of void regions in the flow. -To address these problems, [Adami (2013)](@cite Adami2013) modified the advection velocity and added an extra term to the momentum equation. -The authors introduced the so-called Transport Velocity Formulation (TVF) for WCSPH. [Ramachandran (2019)](@cite Ramachandran2019) applied the TVF -also for the [EDAC](@ref edac) scheme. - -The transport velocity ``\tilde{v}_a`` of particle ``a`` is used to evolve the position of the particle ``r_a`` from one time step to the next by - -```math -\frac{\mathrm{d} r_a}{\mathrm{d}t} = \tilde{v}_a -``` - -and is obtained at every time-step ``\Delta t`` from - -```math -\tilde{v}_a (t + \Delta t) = v_a (t) + \Delta t \left(\frac{\tilde{\mathrm{d}} v_a}{\mathrm{d}t} - \frac{1}{\rho_a} \nabla p_{\text{background}} \right), -``` - -where ``\rho_a`` is the density of particle ``a`` and ``p_{\text{background}}`` is a constant background pressure field. -The tilde in the second term of the right hand side indicates that the material derivative has an advection part. - -The discretized form of the last term is - -```math - -\frac{1}{\rho_a} \nabla p_{\text{background}} \approx -\frac{p_{\text{background}}}{m_a} \sum_b \left(V_a^2 + V_b^2 \right) \nabla_a W_{ab}, -``` - -where ``V_a``, ``V_b`` denote the volume of particles ``a`` and ``b`` respectively. -Note that although in the continuous case ``\nabla p_{\text{background}} = 0``, the discretization is not 0th-order consistent for **non**-uniform particle distribution, -which means that there is a non-vanishing contribution only when particles are disordered. -That also means that ``p_{\text{background}}`` occurs as prefactor to correct the trajectory of a particle resulting in uniform pressure distributions. -Suggested is a background pressure which is in the order of the reference pressure but can be chosen arbitrarily large when the time-step criterion is adjusted. - -The inviscid momentum equation with an additional convection term for a particle moving with ``\tilde{v}`` is - -```math -\frac{\tilde{\mathrm{d}} \left( \rho v \right)}{\mathrm{d}t} = -\nabla p + \nabla \cdot \bm{A}, -``` - - where the tensor ``\bm{A} = \rho v\left(\tilde{v}-v\right)^T`` is a consequence of the modified - advection velocity and can be interpreted as the convection of momentum with the relative velocity ``\tilde{v}-v``. - -The discretized form of the momentum equation for a particle ``a`` reads as - -```math -\frac{\tilde{\mathrm{d}} v_a}{\mathrm{d}t} = \frac{1}{m_a} \sum_b \left(V_a^2 + V_b^2 \right) \left[ -\tilde{p}_{ab} \nabla_a W_{ab} + \frac{1}{2} \left(\bm{A}_a + \bm{A}_b \right) \cdot \nabla_a W_{ab} \right]. -``` - -Here, ``\tilde{p}_{ab}`` is the density-weighted pressure - -```math -\tilde{p}_{ab} = \frac{\rho_b p_a + \rho_a p_b}{\rho_a + \rho_b}, -``` - -with the density ``\rho_a``, ``\rho_b`` and the pressure ``p_a``, ``p_b`` of particles ``a`` and ``b`` respectively. ``\bm{A}_a`` and ``\bm{A}_b`` are the convection tensors for particle ``a`` and ``b`` respectively and are given, e.g. for particle ``a``, as ``\bm{A}_a = \rho v_a\left(\tilde{v}_a-v_a\right)^T``. - -```@autodocs -Modules = [TrixiParticles] -Pages = [joinpath("schemes", "fluid", "transport_velocity.jl")] -``` diff --git a/docs/src/systems/weakly_compressible_sph.md b/docs/src/systems/weakly_compressible_sph.md index 51355961f6..5904ffea2c 100644 --- a/docs/src/systems/weakly_compressible_sph.md +++ b/docs/src/systems/weakly_compressible_sph.md @@ -152,8 +152,74 @@ as explained in [Sun2018](@cite Sun2018) on page 29, right above Equation 9. The ``\delta``-SPH method (WCSPH with density diffusion) together with this formulation of PST is commonly referred to as ``\delta^+``-SPH. -The Particle Shifting Technique can be applied in form -of the [`ParticleShiftingCallback`](@ref). +To apply particle shifting, use the keyword argument `shifting_technique` in the constructor +of a system that supports it. + + +## [Transport Velocity Formulation (TVF)](@id transport_velocity_formulation) + +An alternative formulation is the so-called Transport Velocity Formulation (TVF) +by [Adami (2013)](@cite Adami2013). +[Ramachandran (2019)](@cite Ramachandran2019) applied the TVF also for the [EDAC](@ref edac) +scheme. + +The transport velocity ``\tilde{v}_a`` of particle ``a`` is used to evolve the position +of the particle ``r_a`` from one time step to the next by +```math +\frac{\mathrm{d} r_a}{\mathrm{d}t} = \tilde{v}_a +``` +and is obtained at every time step ``\Delta t`` from +```math +\tilde{v}_a (t + \Delta t) = v_a (t) + \Delta t \left(\frac{\tilde{\mathrm{d}} v_a}{\mathrm{d}t} - \frac{1}{\rho_a} \nabla p_{\text{background}} \right), +``` +where ``\rho_a`` is the density of particle ``a`` and ``p_{\text{background}}`` +is a constant background pressure field. +The tilde in the second term of the right-hand side indicates that the material derivative +has an advection part. + +The discretized form of the last term is +```math + -\frac{1}{\rho_a} \nabla p_{\text{background}} \approx -\frac{p_{\text{background}}}{m_a} \sum_b \left(V_a^2 + V_b^2 \right) \nabla_a W_{ab}, +``` +where ``V_a``, ``V_b`` denote the volume of particles ``a`` and ``b`` respectively. +Note that although in the continuous case ``\nabla p_{\text{background}} = 0``, +the discretization is not 0th-order consistent for **non**-uniform particle distribution, +which means that there is a non-vanishing contribution only when particles are disordered. +That also means that ``p_{\text{background}}`` occurs as pre-factor to correct +the trajectory of a particle resulting in uniform pressure distributions. +Suggested is a background pressure which is in the order of the reference pressure, +but it can be chosen arbitrarily large when the time-step criterion is adjusted. + +The inviscid momentum equation with an additional convection term for a particle +moving with ``\tilde{v}`` is +```math +\frac{\tilde{\mathrm{d}} \left( \rho v \right)}{\mathrm{d}t} = -\nabla p + \nabla \cdot \bm{A}, +``` +where the tensor ``\bm{A} = \rho v\left(\tilde{v}-v\right)^T`` is a consequence +of the modified advection velocity and can be interpreted as the convection of momentum +with the relative velocity ``\tilde{v}-v``. + +The discretized form of the momentum equation for a particle ``a`` reads as +```math +\frac{\tilde{\mathrm{d}} v_a}{\mathrm{d}t} = \frac{1}{m_a} \sum_b \left(V_a^2 + V_b^2 \right) \left[ -\tilde{p}_{ab} \nabla_a W_{ab} + \frac{1}{2} \left(\bm{A}_a + \bm{A}_b \right) \cdot \nabla_a W_{ab} \right]. +``` +Here, ``\tilde{p}_{ab}`` is the density-weighted pressure +```math +\tilde{p}_{ab} = \frac{\rho_b p_a + \rho_a p_b}{\rho_a + \rho_b}, +``` +with the density ``\rho_a``, ``\rho_b`` and the pressure ``p_a``, ``p_b`` of particles ``a`` +and ``b``, respectively. ``\bm{A}_a`` and ``\bm{A}_b`` are the convection tensors +for particle ``a`` and ``b``, respectively, and are given, e.g., for particle ``a``, +as ``\bm{A}_a = \rho v_a\left(\tilde{v}_a-v_a\right)^T``. + +To apply the TVF, use the keyword argument `shifting_technique` in the constructor +of a system that supports it. + +```@autodocs +Modules = [TrixiParticles] +Pages = [joinpath("schemes", "fluid", "shifting_technique.jl")] +``` + ## [Tensile Instability Control](@id tic) From 13d32680ae4086d962c6374f258ab3fd77ef837e Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Sun, 17 Aug 2025 16:33:39 +0300 Subject: [PATCH 07/13] Fix alle example files --- examples/fluid/lid_driven_cavity_2d.jl | 4 ++-- examples/fluid/periodic_array_of_cylinders_2d.jl | 2 +- examples/fluid/periodic_channel_2d.jl | 5 ++--- examples/fluid/pipe_flow_2d.jl | 6 +++--- examples/fluid/taylor_green_vortex_2d.jl | 4 ++-- src/TrixiParticles.jl | 5 ++--- test/examples/examples_fluid.jl | 14 ++++++++++++-- test/systems/edac_system.jl | 2 +- 8 files changed, 25 insertions(+), 17 deletions(-) diff --git a/examples/fluid/lid_driven_cavity_2d.jl b/examples/fluid/lid_driven_cavity_2d.jl index b6582a1701..1b8b15cda4 100644 --- a/examples/fluid/lid_driven_cavity_2d.jl +++ b/examples/fluid/lid_driven_cavity_2d.jl @@ -65,7 +65,7 @@ if wcsph state_equation, smoothing_kernel, pressure_acceleration=TrixiParticles.inter_particle_averaged_pressure, smoothing_length, viscosity=viscosity, - transport_velocity=TransportVelocityAdami(pressure)) + shifting_technique=TransportVelocityAdami(pressure)) else state_equation = nothing density_calculator = ContinuityDensity() @@ -73,7 +73,7 @@ else smoothing_length, density_calculator=density_calculator, sound_speed, viscosity=viscosity, - transport_velocity=TransportVelocityAdami(pressure)) + shifting_technique=TransportVelocityAdami(pressure)) end # ========================================================================================== diff --git a/examples/fluid/periodic_array_of_cylinders_2d.jl b/examples/fluid/periodic_array_of_cylinders_2d.jl index 850c9f2ff9..648a98be11 100644 --- a/examples/fluid/periodic_array_of_cylinders_2d.jl +++ b/examples/fluid/periodic_array_of_cylinders_2d.jl @@ -60,7 +60,7 @@ smoothing_length = 1.2 * particle_spacing smoothing_kernel = SchoenbergQuarticSplineKernel{2}() fluid_system = EntropicallyDampedSPHSystem(fluid, smoothing_kernel, smoothing_length, sound_speed, viscosity=ViscosityAdami(; nu), - transport_velocity=TransportVelocityAdami(pressure), + shifting_technique=TransportVelocityAdami(pressure), acceleration=(acceleration_x, 0.0)) # ========================================================================================== diff --git a/examples/fluid/periodic_channel_2d.jl b/examples/fluid/periodic_channel_2d.jl index 74770c318b..811ba0df2c 100644 --- a/examples/fluid/periodic_channel_2d.jl +++ b/examples/fluid/periodic_channel_2d.jl @@ -48,6 +48,7 @@ viscosity = ArtificialViscosityMonaghan(alpha=0.02, beta=0.0) fluid_system = WeaklyCompressibleSPHSystem(tank.fluid, fluid_density_calculator, state_equation, smoothing_kernel, smoothing_length, viscosity=viscosity, + particle_shifting_technique=nothing, pressure_acceleration=nothing) # ========================================================================================== @@ -76,10 +77,8 @@ ode = semidiscretize(semi, tspan) info_callback = InfoCallback(interval=100) saving_callback = SolutionSavingCallback(dt=0.02, prefix="") -# This can be overwritten with `trixi_include` -extra_callback = nothing -callbacks = CallbackSet(info_callback, saving_callback, extra_callback) +callbacks = CallbackSet(info_callback, saving_callback) # Use a Runge-Kutta method with automatic (error based) time step size control. # Limiting of the maximum stepsize is necessary to prevent crashing. diff --git a/examples/fluid/pipe_flow_2d.jl b/examples/fluid/pipe_flow_2d.jl index 3544ac7a8d..85e3461ccc 100644 --- a/examples/fluid/pipe_flow_2d.jl +++ b/examples/fluid/pipe_flow_2d.jl @@ -74,6 +74,7 @@ viscosity = ViscosityAdami(nu=kinematic_viscosity) fluid_system = EntropicallyDampedSPHSystem(pipe.fluid, smoothing_kernel, smoothing_length, sound_speed, viscosity=viscosity, density_calculator=fluid_density_calculator, + shifting_technique=ParticleShiftingTechnique(), buffer_size=n_buffer_particles) # Alternatively the WCSPH scheme can be used @@ -86,6 +87,7 @@ if wcsph fluid_system = WeaklyCompressibleSPHSystem(pipe.fluid, fluid_density_calculator, state_equation, smoothing_kernel, smoothing_length, viscosity=viscosity, + shifting_technique=ParticleShiftingTechnique(), buffer_size=n_buffer_particles) end @@ -159,12 +161,10 @@ ode = semidiscretize(semi, tspan) info_callback = InfoCallback(interval=100) saving_callback = SolutionSavingCallback(dt=0.02, prefix="") -particle_shifting = ParticleShiftingCallback() extra_callback = nothing -callbacks = CallbackSet(info_callback, saving_callback, UpdateCallback(), - particle_shifting, extra_callback) +callbacks = CallbackSet(info_callback, saving_callback, UpdateCallback(), extra_callback) sol = solve(ode, RDPK3SpFSAL35(), abstol=1e-5, # Default abstol is 1e-6 (may need to be tuned to prevent boundary penetration) diff --git a/examples/fluid/taylor_green_vortex_2d.jl b/examples/fluid/taylor_green_vortex_2d.jl index e6e74fbc4b..8a99c455c6 100644 --- a/examples/fluid/taylor_green_vortex_2d.jl +++ b/examples/fluid/taylor_green_vortex_2d.jl @@ -86,13 +86,13 @@ if wcsph pressure_acceleration=TrixiParticles.inter_particle_averaged_pressure, smoothing_length, viscosity=ViscosityAdami(; nu), - transport_velocity=TransportVelocityAdami(background_pressure)) + shifting_technique=TransportVelocityAdami(background_pressure)) else density_calculator = SummationDensity() fluid_system = EntropicallyDampedSPHSystem(fluid, smoothing_kernel, smoothing_length, sound_speed, density_calculator=density_calculator, - transport_velocity=TransportVelocityAdami(background_pressure), + shifting_technique=TransportVelocityAdami(background_pressure), viscosity=ViscosityAdami(; nu)) end diff --git a/src/TrixiParticles.jl b/src/TrixiParticles.jl index 987e208750..2a53025d46 100644 --- a/src/TrixiParticles.jl +++ b/src/TrixiParticles.jl @@ -65,10 +65,9 @@ export WeaklyCompressibleSPHSystem, EntropicallyDampedSPHSystem, TotalLagrangian BoundarySPHSystem, DEMSystem, BoundaryDEMSystem, OpenBoundarySPHSystem export BoundaryZone, InFlow, OutFlow, BidirectionalFlow export InfoCallback, SolutionSavingCallback, DensityReinitializationCallback, - PostprocessCallback, StepsizeCallback, UpdateCallback, SteadyStateReachedCallback, - ParticleShiftingCallback + PostprocessCallback, StepsizeCallback, UpdateCallback, SteadyStateReachedCallback export ContinuityDensity, SummationDensity -export PenaltyForceGanzenmueller, TransportVelocityAdami +export PenaltyForceGanzenmueller, TransportVelocityAdami, ParticleShiftingTechnique export SchoenbergCubicSplineKernel, SchoenbergQuarticSplineKernel, SchoenbergQuinticSplineKernel, GaussianKernel, WendlandC2Kernel, WendlandC4Kernel, WendlandC6Kernel, SpikyKernel, Poly6Kernel diff --git a/test/examples/examples_fluid.jl b/test/examples/examples_fluid.jl index ede9e6bbf9..6432f4841d 100644 --- a/test/examples/examples_fluid.jl +++ b/test/examples/examples_fluid.jl @@ -273,7 +273,17 @@ joinpath(examples_dir(), "fluid", "periodic_channel_2d.jl"), tspan=(0.0, 0.2), - extra_callback=ParticleShiftingCallback()) + shifting_technique=ParticleShiftingTechnique()) + @test sol.retcode == ReturnCode.Success + @test count_rhs_allocations(sol, semi) == 0 + end + + @trixi_testset "fluid/periodic_channel_2d.jl with TVF" begin + @trixi_test_nowarn trixi_include(@__MODULE__, + joinpath(examples_dir(), "fluid", + "periodic_channel_2d.jl"), + tspan=(0.0, 0.2), + shifting_technique=TransportVelocityAdami(50_000.0)) @test sol.retcode == ReturnCode.Success @test count_rhs_allocations(sol, semi) == 0 end @@ -283,7 +293,7 @@ joinpath(examples_dir(), "fluid", "periodic_channel_2d.jl"), tspan=(0.0, 0.2), - extra_callback=ParticleShiftingCallback(), + shifting_technique=ParticleShiftingTechnique(), pressure_acceleration=tensile_instability_control) @test sol.retcode == ReturnCode.Success @test count_rhs_allocations(sol, semi) == 0 diff --git a/test/systems/edac_system.jl b/test/systems/edac_system.jl index cdebb91768..1d31afdcff 100644 --- a/test/systems/edac_system.jl +++ b/test/systems/edac_system.jl @@ -225,7 +225,7 @@ names = ["No TVF", "TransportVelocityAdami"] @testset "$(names[i])" for i in eachindex(transport_velocity) system = EntropicallyDampedSPHSystem(fluid, smoothing_kernel, - transport_velocity=transport_velocity[i], + shifting_technique=transport_velocity[i], average_pressure_reduction=true, smoothing_length, 0.0) semi = Semidiscretization(system) From 2a0ae285284805afb5102d530090ff67535dba84 Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Sun, 17 Aug 2025 16:42:30 +0300 Subject: [PATCH 08/13] Fix tests --- src/general/semidiscretization.jl | 2 +- test/systems/edac_system.jl | 2 +- test/systems/wcsph_system.jl | 4 ++-- 3 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/general/semidiscretization.jl b/src/general/semidiscretization.jl index a9c78b3349..9532d812ab 100644 --- a/src/general/semidiscretization.jl +++ b/src/general/semidiscretization.jl @@ -485,7 +485,7 @@ end @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 + # This is zero unless a shifting technique is used delta_v_ = delta_v(system, particle) for i in 1:ndims(system) diff --git a/test/systems/edac_system.jl b/test/systems/edac_system.jl index 1d31afdcff..ed37b9ea55 100644 --- a/test/systems/edac_system.jl +++ b/test/systems/edac_system.jl @@ -140,7 +140,7 @@ │ viscosity: …………………………………………………… Nothing │ │ ν₍EDAC₎: ………………………………………………………… ≈ 0.226 │ │ smoothing kernel: ………………………………… Val │ - │ tansport velocity formulation: Nothing │ + │ shifting technique: …………………………… nothing │ │ average pressure reduction: ……… no │ │ acceleration: …………………………………………… [0.0, 0.0] │ │ surface tension: …………………………………… nothing │ diff --git a/test/systems/wcsph_system.jl b/test/systems/wcsph_system.jl index cf6020d068..99f7152db6 100644 --- a/test/systems/wcsph_system.jl +++ b/test/systems/wcsph_system.jl @@ -199,7 +199,7 @@ smoothing_length, density_diffusion=density_diffusion) - show_compact = "WeaklyCompressibleSPHSystem{2}(SummationDensity(), nothing, Val{:state_equation}(), Val{:smoothing_kernel}(), nothing, Val{:density_diffusion}(), nothing, nothing, [0.0, 0.0], nothing) with 2 particles" + show_compact = "WeaklyCompressibleSPHSystem{2}(SummationDensity(), nothing, Val{:state_equation}(), Val{:smoothing_kernel}(), nothing, Val{:density_diffusion}(), nothing, nothing, nothing, [0.0, 0.0], nothing) with 2 particles" @test repr(system) == show_compact show_box = """ ┌──────────────────────────────────────────────────────────────────────────────────────────────────┐ @@ -211,8 +211,8 @@ │ state equation: ……………………………………… Val │ │ smoothing kernel: ………………………………… Val │ │ viscosity: …………………………………………………… nothing │ - │ tansport velocity formulation: Nothing │ │ density diffusion: ……………………………… Val{:density_diffusion}() │ + │ shifting technique: …………………………… nothing │ │ surface tension: …………………………………… nothing │ │ surface normal method: …………………… nothing │ │ acceleration: …………………………………………… [0.0, 0.0] │ From b5d1492bbce32768e0eb60027bb867c4a539139b Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Mon, 18 Aug 2025 12:10:59 +0300 Subject: [PATCH 09/13] Fix periodic channel --- examples/fluid/periodic_channel_2d.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/fluid/periodic_channel_2d.jl b/examples/fluid/periodic_channel_2d.jl index 811ba0df2c..8d6ef4458f 100644 --- a/examples/fluid/periodic_channel_2d.jl +++ b/examples/fluid/periodic_channel_2d.jl @@ -28,7 +28,7 @@ initial_fluid_size = tank_size initial_velocity = (1.0, 0.0) fluid_density = 1000.0 -sound_speed = initial_velocity[1] +sound_speed = 10 * initial_velocity[1] state_equation = StateEquationCole(; sound_speed, reference_density=fluid_density, exponent=7) @@ -48,7 +48,7 @@ viscosity = ArtificialViscosityMonaghan(alpha=0.02, beta=0.0) fluid_system = WeaklyCompressibleSPHSystem(tank.fluid, fluid_density_calculator, state_equation, smoothing_kernel, smoothing_length, viscosity=viscosity, - particle_shifting_technique=nothing, + shifting_technique=nothing, pressure_acceleration=nothing) # ========================================================================================== From 77e20e1148c740cdf8132dcc4465ac142f3cb639 Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Mon, 18 Aug 2025 12:14:47 +0300 Subject: [PATCH 10/13] Fix docs --- docs/src/systems/weakly_compressible_sph.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/systems/weakly_compressible_sph.md b/docs/src/systems/weakly_compressible_sph.md index 5904ffea2c..62a1bc8d0c 100644 --- a/docs/src/systems/weakly_compressible_sph.md +++ b/docs/src/systems/weakly_compressible_sph.md @@ -217,7 +217,7 @@ of a system that supports it. ```@autodocs Modules = [TrixiParticles] -Pages = [joinpath("schemes", "fluid", "shifting_technique.jl")] +Pages = [joinpath("schemes", "fluid", "shifting_techniques.jl")] ``` From 0209d65ca390bfa9426c08221a2a5735f4f381ce Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Mon, 18 Aug 2025 12:27:46 +0300 Subject: [PATCH 11/13] Update news --- NEWS.md | 25 ++++++++++++++++++------ src/schemes/fluid/shifting_techniques.jl | 5 +++++ 2 files changed, 24 insertions(+), 6 deletions(-) diff --git a/NEWS.md b/NEWS.md index f1940cbc46..d4e2ccdd7d 100644 --- a/NEWS.md +++ b/NEWS.md @@ -4,14 +4,27 @@ TrixiParticles.jl follows the interpretation of [semantic versioning (semver)](https://julialang.github.io/Pkg.jl/dev/compatibility/#Version-specifier-format-1) used in the Julia ecosystem. Notable changes will be documented in this file for human readability. +## Version 0.4 + +### API Changes + +- Combined transport velocity formulation (TVF) and particle shifting technique (PST) into + one unified framework. + The keyword argument `transport_velocity` now changed to `shifting_technique`. + The `ParticleShiftingCallback` has been removed. To use PST, use the `UpdateCallback` + instead, and pass `shifting_technique=ParticleShiftingTechnique()` to the system. + +- Renamed the keyword argument `tlsph` to `place_on_shell` for `ParticlePackingSystem`, + `sample_boundary`, `extrude_geometry`, `RectangularShape`, and `SphereShape`. + ## Version 0.3.1 ### Features -- **Simplified SGS Viscosity Models**: Added ViscosityMorrisSGS and ViscosityAdamiSGS, +- **Simplified SGS Viscosity Models**: Added ViscosityMorrisSGS and ViscosityAdamiSGS, which implement a simplified Smagorinsky-type sub-grid-scale viscosity. (#753) -- **Multithreaded Integration Array**: Introduced a new array type for CPU backends +- **Multithreaded Integration Array**: Introduced a new array type for CPU backends that enables multithreaded broadcasting, delivering speed-ups of up to 5× on systems with many threads when combined with thread pinning. (#722) @@ -21,17 +34,17 @@ used in the Julia ecosystem. Notable changes will be documented in this file for - **DXF file format support**: Import complex geometries using the DXF file format. (#821) - **Improved Plane interpolation**: Massively improved interpolation performance for planes (#763). - + ### GPU - Make PST GPU-compatible (#813). - + - Make open boundaries GPU-compatible (#773). - + - Make interpolation GPU-compatible (#812). ### Important Bugfixes - + - Fix validation setups (#801). - Calculate interpolated density instead of computed density when using interpolation (#808). diff --git a/src/schemes/fluid/shifting_techniques.jl b/src/schemes/fluid/shifting_techniques.jl index 13c78a4f9c..98775b3be0 100644 --- a/src/schemes/fluid/shifting_techniques.jl +++ b/src/schemes/fluid/shifting_techniques.jl @@ -174,6 +174,11 @@ 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. + +!!! warning + The Transport Velocity Formulation needs to be disabled close to the free surface + and therefore requires a free surface detection method. This is not yet implemented. + **This technique cannot be used in a free surface simulation.** """ struct TransportVelocityAdami{T <: Real} <: AbstractShiftingTechnique background_pressure::T From e8146dcfd900b65a91986ab2efbf6931d8c58c04 Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Thu, 21 Aug 2025 12:03:59 +0300 Subject: [PATCH 12/13] Fix tests --- test/examples/examples_fluid.jl | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/test/examples/examples_fluid.jl b/test/examples/examples_fluid.jl index 6432f4841d..a84a69317d 100644 --- a/test/examples/examples_fluid.jl +++ b/test/examples/examples_fluid.jl @@ -273,7 +273,8 @@ joinpath(examples_dir(), "fluid", "periodic_channel_2d.jl"), tspan=(0.0, 0.2), - shifting_technique=ParticleShiftingTechnique()) + shifting_technique=ParticleShiftingTechnique(), + extra_callback=UpdateCallback()) @test sol.retcode == ReturnCode.Success @test count_rhs_allocations(sol, semi) == 0 end @@ -283,7 +284,8 @@ joinpath(examples_dir(), "fluid", "periodic_channel_2d.jl"), tspan=(0.0, 0.2), - shifting_technique=TransportVelocityAdami(50_000.0)) + shifting_technique=TransportVelocityAdami(50_000.0), + extra_callback=UpdateCallback()) @test sol.retcode == ReturnCode.Success @test count_rhs_allocations(sol, semi) == 0 end @@ -294,7 +296,8 @@ "periodic_channel_2d.jl"), tspan=(0.0, 0.2), shifting_technique=ParticleShiftingTechnique(), - pressure_acceleration=tensile_instability_control) + pressure_acceleration=tensile_instability_control, + extra_callback=UpdateCallback()) @test sol.retcode == ReturnCode.Success @test count_rhs_allocations(sol, semi) == 0 end From 844145e778e6894efa41586e9e2023826f38a194 Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Thu, 21 Aug 2025 12:05:37 +0300 Subject: [PATCH 13/13] Fix example file --- examples/fluid/periodic_channel_2d.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/examples/fluid/periodic_channel_2d.jl b/examples/fluid/periodic_channel_2d.jl index 8d6ef4458f..9cd4bd95a5 100644 --- a/examples/fluid/periodic_channel_2d.jl +++ b/examples/fluid/periodic_channel_2d.jl @@ -77,8 +77,10 @@ ode = semidiscretize(semi, tspan) info_callback = InfoCallback(interval=100) saving_callback = SolutionSavingCallback(dt=0.02, prefix="") +# This can be overwritten with `trixi_include` +extra_callback = nothing -callbacks = CallbackSet(info_callback, saving_callback) +callbacks = CallbackSet(info_callback, saving_callback, extra_callback) # Use a Runge-Kutta method with automatic (error based) time step size control. # Limiting of the maximum stepsize is necessary to prevent crashing.