From a5c158c0e10841cb84a7e907edf02bdae08a9c92 Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Mon, 11 May 2026 17:41:02 +0200 Subject: [PATCH 1/5] Fix density reinitialization callback --- src/callbacks/density_reinit.jl | 30 ++++++++++++++--- test/callbacks/callbacks.jl | 1 + test/callbacks/density_reinit.jl | 57 ++++++++++++++++++++++++++++++++ 3 files changed, 83 insertions(+), 5 deletions(-) create mode 100644 test/callbacks/density_reinit.jl diff --git a/src/callbacks/density_reinit.jl b/src/callbacks/density_reinit.jl index ded47149f9..302ae6177b 100644 --- a/src/callbacks/density_reinit.jl +++ b/src/callbacks/density_reinit.jl @@ -9,8 +9,9 @@ Callback to reinitialize the density field when using [`ContinuityDensity`](@ref of integration time. - `reinit_initial_solution`: Reinitialize the initial solution (default=false) """ -mutable struct DensityReinitializationCallback{I} +mutable struct DensityReinitializationCallback{I, PS} interval::I + particle_system::PS last_t::Float64 reinit_initial_solution::Bool end @@ -29,7 +30,7 @@ function Base.show(io::IO, ::MIME"text/plain", show(io, cb) else callback = cb.affect! - setup = [ + setup = Pair{String, Any}[ "interval" => callback.interval, "reinit_initial_solution" => callback.reinit_initial_solution ] @@ -53,7 +54,8 @@ function DensityReinitializationCallback(particle_system; interval::Integer=0, d last_t = -Inf - reinit_cb = DensityReinitializationCallback(interval, last_t, reinit_initial_solution) + reinit_cb = DensityReinitializationCallback(interval, particle_system, last_t, + reinit_initial_solution) return DiscreteCallback(reinit_cb, reinit_cb, save_positions=(false, false), initialize=(initialize_reinit_cb!)) @@ -91,7 +93,7 @@ end function (reinit_callback::DensityReinitializationCallback)(u, t, integrator) (; interval, last_t) = reinit_callback - return (t - last_t) > interval + return (t - last_t) >= interval end # affect! @@ -99,7 +101,25 @@ function (reinit_callback::DensityReinitializationCallback)(integrator) vu_ode = integrator.u semi = integrator.p - @trixi_timeit timer() "reinit density" reinit_density!(vu_ode, semi) + @trixi_timeit timer() "reinit density" reinitialize_density!(reinit_callback, vu_ode, + semi) reinit_callback.last_t = integrator.t + + # Tell OrdinaryDiffEq that `integrator.u` has been modified + u_modified!(integrator, true) + + return integrator +end + +function reinitialize_density!(reinit_callback::DensityReinitializationCallback, + vu_ode, semi) + (; particle_system) = reinit_callback + v_ode, u_ode = vu_ode.x + v = wrap_v(v_ode, particle_system, semi) + u = wrap_u(u_ode, particle_system, semi) + + reinit_density!(particle_system, v, u, v_ode, u_ode, semi) + + return reinit_callback end diff --git a/test/callbacks/callbacks.jl b/test/callbacks/callbacks.jl index b54fd19c3f..e729a3f0e3 100644 --- a/test/callbacks/callbacks.jl +++ b/test/callbacks/callbacks.jl @@ -55,5 +55,6 @@ include("postprocess.jl") include("update.jl") include("solution_saving.jl") + include("density_reinit.jl") include("steady_state_reached.jl") end diff --git a/test/callbacks/density_reinit.jl b/test/callbacks/density_reinit.jl new file mode 100644 index 0000000000..1c4ec555fb --- /dev/null +++ b/test/callbacks/density_reinit.jl @@ -0,0 +1,57 @@ +@testset verbose=true "DensityReinitializationCallback" begin + struct MockDensityReinitSystem + density_calculator + name::Symbol + end + + density_reinit_calls = Symbol[] + + TrixiParticles.wrap_v(v_ode, system::MockDensityReinitSystem, semi) = (:v, system.name) + TrixiParticles.wrap_u(u_ode, system::MockDensityReinitSystem, semi) = (:u, system.name) + + function TrixiParticles.reinit_density!(system::MockDensityReinitSystem, v, u, + v_ode, u_ode, semi) + push!(density_reinit_calls, system.name) + return system + end + + @testset verbose=true "show" begin + system = MockDensityReinitSystem(nothing, :fluid) + callback = DensityReinitializationCallback(system; interval=10) + + show_compact = "DensityReinitializationCallback(interval=10, reinit_initial_solution=true)" + @test repr(callback) == show_compact + + show_box = """ + ┌──────────────────────────────────────────────────────────────────────────────────────────────────┐ + │ DensityReinitializationCallback │ + │ ═══════════════════════════════ │ + │ interval: ……………………………………………………… 10 │ + │ reinit_initial_solution: ……………… true │ + └──────────────────────────────────────────────────────────────────────────────────────────────────┘""" + @test repr("text/plain", callback) == show_box + end + + @testset verbose=true "dt condition" begin + system = MockDensityReinitSystem(nothing, :fluid) + callback = DensityReinitializationCallback(system; dt=0.1).affect! + callback.last_t = 0.0 + + @test !callback(nothing, 0.099, nothing) + @test callback(nothing, 0.1, nothing) + end + + @testset verbose=true "selected system" begin + empty!(density_reinit_calls) + + system1 = MockDensityReinitSystem(nothing, :fluid1) + system2 = MockDensityReinitSystem(nothing, :fluid2) + callback = DensityReinitializationCallback(system1; interval=1).affect! + vu_ode = (; x=(:v_ode, :u_ode)) + semi = (; systems=(system1, system2)) + + TrixiParticles.reinitialize_density!(callback, vu_ode, semi) + + @test density_reinit_calls == [:fluid1] + end +end From a01dbc840991d09325667efadb5c3b5fe8f7690c Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Tue, 12 May 2026 13:11:37 +0200 Subject: [PATCH 2/5] format --- src/callbacks/density_reinit.jl | 6 ++---- test/callbacks/density_reinit.jl | 2 +- 2 files changed, 3 insertions(+), 5 deletions(-) diff --git a/src/callbacks/density_reinit.jl b/src/callbacks/density_reinit.jl index 302ae6177b..338ca1bcb9 100644 --- a/src/callbacks/density_reinit.jl +++ b/src/callbacks/density_reinit.jl @@ -30,10 +30,8 @@ function Base.show(io::IO, ::MIME"text/plain", show(io, cb) else callback = cb.affect! - setup = Pair{String, Any}[ - "interval" => callback.interval, - "reinit_initial_solution" => callback.reinit_initial_solution - ] + setup = Pair{String, Any}["interval" => callback.interval, + "reinit_initial_solution" => callback.reinit_initial_solution] summary_box(io, "DensityReinitializationCallback", setup) end end diff --git a/test/callbacks/density_reinit.jl b/test/callbacks/density_reinit.jl index 1c4ec555fb..809fa53bbb 100644 --- a/test/callbacks/density_reinit.jl +++ b/test/callbacks/density_reinit.jl @@ -1,6 +1,6 @@ @testset verbose=true "DensityReinitializationCallback" begin struct MockDensityReinitSystem - density_calculator + density_calculator::Any name::Symbol end From 2eaa387e61bb1d06099db6a9d609bfd141dcd7bb Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Sun, 17 May 2026 21:42:40 +0200 Subject: [PATCH 3/5] fix --- examples/fluid/dam_break_2d.jl | 2 +- src/callbacks/density_reinit.jl | 168 +++++++++++++++++++++++++++++-- test/callbacks/density_reinit.jl | 64 ++++++++++++ 3 files changed, 222 insertions(+), 12 deletions(-) diff --git a/examples/fluid/dam_break_2d.jl b/examples/fluid/dam_break_2d.jl index c14119cb60..58b40c78be 100644 --- a/examples/fluid/dam_break_2d.jl +++ b/examples/fluid/dam_break_2d.jl @@ -112,7 +112,7 @@ extra_callback2 = nothing use_reinit = false density_reinit_cb = use_reinit ? - DensityReinitializationCallback(semi.systems[1], interval=10) : + DensityReinitializationCallback(system_index=1, interval=10) : nothing stepsize_callback = StepsizeCallback(cfl=0.9) diff --git a/src/callbacks/density_reinit.jl b/src/callbacks/density_reinit.jl index 2e02f44996..2211f829c3 100644 --- a/src/callbacks/density_reinit.jl +++ b/src/callbacks/density_reinit.jl @@ -1,17 +1,25 @@ """ - DensityReinitializationCallback(; interval::Integer=0, dt=0.0) + DensityReinitializationCallback(particle_system; interval::Integer=0, dt=0.0) + DensityReinitializationCallback(system_index::Integer; interval::Integer=0, dt=0.0) + DensityReinitializationCallback(system_indices; interval::Integer=0, dt=0.0) + DensityReinitializationCallback(; system_index::Integer, interval::Integer=0, dt=0.0) + DensityReinitializationCallback(; system_indices, interval::Integer=0, dt=0.0) Callback to reinitialize the density field when using [`ContinuityDensity`](@ref) [Panizzo2007](@cite). +Pass `system_index` or `system_indices` to select systems by their position in the +semidiscretization. This is robust when [`semidiscretize`](@ref) replaces systems internally. + # Keywords - `interval=0`: Reinitialize the density every `interval` time steps. - `dt`: Reinitialize the density in regular intervals of `dt` in terms of integration time. - `reinit_initial_solution`: Reinitialize the initial solution (default=false) """ -mutable struct DensityReinitializationCallback{I, PS} +mutable struct DensityReinitializationCallback{I, PS, SI} interval::I particle_system::PS + system_indices::SI last_t::Float64 reinit_initial_solution::Bool end @@ -46,28 +54,90 @@ function DensityReinitializationCallback(particle_system; interval::Integer=0, d interval = Float64(dt) end - if particle_system.density_calculator isa SummationDensity - throw(ArgumentError("density reinitialization doesn't provide any advantage for summation density")) - end + check_density_reinit_system(particle_system) last_t = -Inf - reinit_cb = DensityReinitializationCallback(interval, particle_system, last_t, + reinit_cb = DensityReinitializationCallback(interval, particle_system, nothing, last_t, reinit_initial_solution) return DiscreteCallback(reinit_cb, reinit_cb, save_positions=(false, false), initialize=(initialize_reinit_cb!)) end +function DensityReinitializationCallback(system_index::Integer; interval::Integer=0, dt=0.0, + reinit_initial_solution=true) + return reinit_callback_from_indices((system_index,); interval, dt, + reinit_initial_solution) +end + +function DensityReinitializationCallback(system_indices::Tuple{Vararg{Integer}}; + interval::Integer=0, dt=0.0, + reinit_initial_solution=true) + return reinit_callback_from_indices(system_indices; interval, dt, reinit_initial_solution) +end + +function DensityReinitializationCallback(system_indices::AbstractVector{<:Integer}; + interval::Integer=0, dt=0.0, + reinit_initial_solution=true) + return reinit_callback_from_indices(system_indices; interval, dt, reinit_initial_solution) +end + +function DensityReinitializationCallback(system_indices::AbstractRange{<:Integer}; + interval::Integer=0, dt=0.0, + reinit_initial_solution=true) + return reinit_callback_from_indices(system_indices; interval, dt, reinit_initial_solution) +end + +function reinit_callback_from_indices(system_indices; interval::Integer=0, dt=0.0, + reinit_initial_solution=true) + if dt > 0 && interval > 0 + error("Setting both interval and dt is not supported!") + end + + if dt > 0 + interval = Float64(dt) + end + + last_t = -Inf + system_indices_ = normalize_reinit_system_indices(system_indices) + + reinit_cb = DensityReinitializationCallback(interval, nothing, system_indices_, + last_t, reinit_initial_solution) + + return DiscreteCallback(reinit_cb, reinit_cb, save_positions=(false, false), + initialize=(initialize_reinit_cb!)) +end + +function DensityReinitializationCallback(; system_index=nothing, system_indices=nothing, + interval::Integer=0, + dt=0.0, reinit_initial_solution=true) + if !isnothing(system_index) && !isnothing(system_indices) + throw(ArgumentError("pass either `system_index` or `system_indices`, not both")) + end + + if isnothing(system_index) && isnothing(system_indices) + throw(ArgumentError("pass `system_index` or `system_indices`")) + end + + indices = isnothing(system_index) ? system_indices : system_index + + return reinit_callback_from_indices(indices; interval, dt, reinit_initial_solution) +end + function initialize_reinit_cb!(cb, u, t, integrator) initialize_reinit_cb!(cb.affect!, u, t, integrator) end function initialize_reinit_cb!(cb::DensityReinitializationCallback, u, t, integrator) + semi = integrator.p.semi + foreach_reinit_system(cb, semi) do system + check_density_reinit_system(system) + end + # Reinitialize initial solution if cb.reinit_initial_solution # Update systems to compute quantities like density and pressure. - semi = integrator.p.semi v_ode, u_ode = u.x update_systems_and_nhs(v_ode, u_ode, semi, t) @@ -112,12 +182,88 @@ end function reinitialize_density!(reinit_callback::DensityReinitializationCallback, vu_ode, semi) - (; particle_system) = reinit_callback v_ode, u_ode = vu_ode.x - v = wrap_v(v_ode, particle_system, semi) - u = wrap_u(u_ode, particle_system, semi) - reinit_density!(particle_system, v, u, v_ode, u_ode, semi) + foreach_reinit_system(reinit_callback, semi) do particle_system + check_density_reinit_system(particle_system) + v = wrap_v(v_ode, particle_system, semi) + u = wrap_u(u_ode, particle_system, semi) + + reinit_density!(particle_system, v, u, v_ode, u_ode, semi) + end return reinit_callback end + +function foreach_reinit_system(f, reinit_callback::DensityReinitializationCallback, semi) + (; particle_system, system_indices) = reinit_callback + + if !isnothing(system_indices) + foreach(system_indices) do system_index + f(current_reinit_system(system_index, semi)) + end + + return nothing + end + + system_index = findfirst(system -> system === particle_system, semi.systems) + + if isnothing(system_index) + throw(ArgumentError("system is not in the semidiscretization. This can happen " * + "when `semidiscretize` replaced the system. Pass a system " * + "index to `DensityReinitializationCallback` to select the " * + "system independently of replacement.")) + end + + f(semi.systems[system_index]) + + return nothing +end + +function current_reinit_system(system_index, semi) + if system_index > length(semi.systems) + throw(ArgumentError("system index $system_index is out of bounds for a " * + "semidiscretization with $(length(semi.systems)) systems")) + end + + return semi.systems[system_index] +end + +function check_density_reinit_system(particle_system) + if !hasproperty(particle_system, :density_calculator) + throw(ArgumentError("density reinitialization requires a system with a density calculator")) + end + + if particle_system.density_calculator isa SummationDensity + throw(ArgumentError("density reinitialization doesn't provide any advantage for summation density")) + end + + return particle_system +end + +function normalize_reinit_system_indices(system_index::Integer) + return (normalize_reinit_system_index(system_index),) +end + +function normalize_reinit_system_indices(system_indices) + system_indices_ = Tuple(normalize_reinit_system_index(system_index) + for system_index in system_indices) + + if isempty(system_indices_) + throw(ArgumentError("`system_indices` must not be empty")) + end + + if !allunique(system_indices_) + throw(ArgumentError("`system_indices` must be unique")) + end + + return system_indices_ +end + +function normalize_reinit_system_index(system_index::Integer) + if system_index < 1 + throw(ArgumentError("system indices must be positive")) + end + + return Int(system_index) +end diff --git a/test/callbacks/density_reinit.jl b/test/callbacks/density_reinit.jl index 809fa53bbb..1587378569 100644 --- a/test/callbacks/density_reinit.jl +++ b/test/callbacks/density_reinit.jl @@ -4,6 +4,10 @@ name::Symbol end + struct MockNoDensityReinitSystem + name::Symbol + end + density_reinit_calls = Symbol[] TrixiParticles.wrap_v(v_ode, system::MockDensityReinitSystem, semi) = (:v, system.name) @@ -54,4 +58,64 @@ @test density_reinit_calls == [:fluid1] end + + @testset verbose=true "selected system after semidiscretization replacement" begin + empty!(density_reinit_calls) + + replacement_system1 = MockDensityReinitSystem(Ref(:replacement1), :replacement1) + replacement_system2 = MockDensityReinitSystem(Ref(:replacement2), :replacement2) + + callback = DensityReinitializationCallback(system_index=1, interval=1).affect! + vu_ode = (; x=(:v_ode, :u_ode)) + semi_replaced = (; systems=(replacement_system1, replacement_system2)) + + TrixiParticles.reinitialize_density!(callback, vu_ode, semi_replaced) + + @test density_reinit_calls == [:replacement1] + end + + @testset verbose=true "selected systems after semidiscretization replacement" begin + empty!(density_reinit_calls) + + replacement_system1 = MockDensityReinitSystem(Ref(:replacement1), :replacement1) + replacement_system2 = MockDensityReinitSystem(Ref(:replacement2), :replacement2) + replacement_system3 = MockDensityReinitSystem(Ref(:replacement3), :replacement3) + + callback = DensityReinitializationCallback(system_indices=(3, 1), + interval=1).affect! + vu_ode = (; x=(:v_ode, :u_ode)) + semi_replaced = (; systems=(replacement_system1, replacement_system2, + replacement_system3)) + + TrixiParticles.reinitialize_density!(callback, vu_ode, semi_replaced) + + @test density_reinit_calls == [:replacement3, :replacement1] + end + + @testset verbose=true "system index validation" begin + @test_throws ArgumentError DensityReinitializationCallback(system_index=0, + interval=1) + @test_throws ArgumentError DensityReinitializationCallback(system_indices=(), + interval=1) + @test_throws ArgumentError DensityReinitializationCallback(system_indices=(1, 1), + interval=1) + + callback = DensityReinitializationCallback(system_index=3, interval=1).affect! + vu_ode = (; x=(:v_ode, :u_ode)) + semi = (; systems=(MockDensityReinitSystem(nothing, :fluid),)) + + @test_throws ArgumentError TrixiParticles.reinitialize_density!(callback, vu_ode, + semi) + + callback = DensityReinitializationCallback(system_index=1, interval=1).affect! + semi = (; systems=(MockDensityReinitSystem(SummationDensity(), :fluid),)) + + @test_throws ArgumentError TrixiParticles.reinitialize_density!(callback, vu_ode, + semi) + + semi = (; systems=(MockNoDensityReinitSystem(:boundary),)) + + @test_throws ArgumentError TrixiParticles.reinitialize_density!(callback, vu_ode, + semi) + end end From 5ae72f322b927212b716fd2f778672e951de66ba Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Mon, 18 May 2026 12:38:58 +0200 Subject: [PATCH 4/5] fix --- examples/fluid/dam_break_2d.jl | 3 +- src/callbacks/density_reinit.jl | 179 +++++++------------------------ test/callbacks/density_reinit.jl | 106 ++++++++++-------- 3 files changed, 102 insertions(+), 186 deletions(-) diff --git a/examples/fluid/dam_break_2d.jl b/examples/fluid/dam_break_2d.jl index 58b40c78be..ee6f837dfa 100644 --- a/examples/fluid/dam_break_2d.jl +++ b/examples/fluid/dam_break_2d.jl @@ -111,8 +111,9 @@ extra_callback = nothing extra_callback2 = nothing use_reinit = false +# Use the system stored in `ode.p.semi`, since `semidiscretize` can replace systems. density_reinit_cb = use_reinit ? - DensityReinitializationCallback(system_index=1, interval=10) : + DensityReinitializationCallback(ode.p.semi.systems[1]; interval=10) : nothing stepsize_callback = StepsizeCallback(cfl=0.9) diff --git a/src/callbacks/density_reinit.jl b/src/callbacks/density_reinit.jl index 2211f829c3..de9a5bbb01 100644 --- a/src/callbacks/density_reinit.jl +++ b/src/callbacks/density_reinit.jl @@ -1,30 +1,12 @@ -""" - DensityReinitializationCallback(particle_system; interval::Integer=0, dt=0.0) - DensityReinitializationCallback(system_index::Integer; interval::Integer=0, dt=0.0) - DensityReinitializationCallback(system_indices; interval::Integer=0, dt=0.0) - DensityReinitializationCallback(; system_index::Integer, interval::Integer=0, dt=0.0) - DensityReinitializationCallback(; system_indices, interval::Integer=0, dt=0.0) - -Callback to reinitialize the density field when using [`ContinuityDensity`](@ref) [Panizzo2007](@cite). - -Pass `system_index` or `system_indices` to select systems by their position in the -semidiscretization. This is robust when [`semidiscretize`](@ref) replaces systems internally. - -# Keywords -- `interval=0`: Reinitialize the density every `interval` time steps. -- `dt`: Reinitialize the density in regular intervals of `dt` in terms - of integration time. -- `reinit_initial_solution`: Reinitialize the initial solution (default=false) -""" -mutable struct DensityReinitializationCallback{I, PS, SI} +mutable struct DensityReinitializationCallbackAffect{I, PS} interval::I - particle_system::PS - system_indices::SI + system::PS last_t::Float64 reinit_initial_solution::Bool end -function Base.show(io::IO, cb::DiscreteCallback{<:Any, <:DensityReinitializationCallback}) +function Base.show(io::IO, + cb::DiscreteCallback{<:Any, <:DensityReinitializationCallbackAffect}) @nospecialize cb # reduce precompilation time callback = cb.affect! print(io, "DensityReinitializationCallback(interval=", callback.interval, @@ -32,7 +14,7 @@ function Base.show(io::IO, cb::DiscreteCallback{<:Any, <:DensityReinitialization end function Base.show(io::IO, ::MIME"text/plain", - cb::DiscreteCallback{<:Any, <:DensityReinitializationCallback}) + cb::DiscreteCallback{<:Any, <:DensityReinitializationCallbackAffect}) @nospecialize cb # reduce precompilation time if get(io, :compact, false) show(io, cb) @@ -44,53 +26,24 @@ function Base.show(io::IO, ::MIME"text/plain", end end -function DensityReinitializationCallback(particle_system; interval::Integer=0, dt=0.0, - reinit_initial_solution=true) - if dt > 0 && interval > 0 - error("Setting both interval and dt is not supported!") - end - - if dt > 0 - interval = Float64(dt) - end - - check_density_reinit_system(particle_system) - - last_t = -Inf - - reinit_cb = DensityReinitializationCallback(interval, particle_system, nothing, last_t, - reinit_initial_solution) - - return DiscreteCallback(reinit_cb, reinit_cb, save_positions=(false, false), - initialize=(initialize_reinit_cb!)) -end - -function DensityReinitializationCallback(system_index::Integer; interval::Integer=0, dt=0.0, - reinit_initial_solution=true) - return reinit_callback_from_indices((system_index,); interval, dt, - reinit_initial_solution) -end +""" + DensityReinitializationCallback(system; interval::Integer=0, dt=0.0, + reinit_initial_solution=true) -function DensityReinitializationCallback(system_indices::Tuple{Vararg{Integer}}; - interval::Integer=0, dt=0.0, - reinit_initial_solution=true) - return reinit_callback_from_indices(system_indices; interval, dt, reinit_initial_solution) -end +Callback to reinitialize the density field when using [`ContinuityDensity`](@ref) [Panizzo2007](@cite). -function DensityReinitializationCallback(system_indices::AbstractVector{<:Integer}; - interval::Integer=0, dt=0.0, - reinit_initial_solution=true) - return reinit_callback_from_indices(system_indices; interval, dt, reinit_initial_solution) -end +Pass a system to reinitialize its density field. The system must be the system stored +in the semidiscretization used by the integrator. If [`semidiscretize`](@ref) creates +a copy of the system, pass the corresponding system from `ode.p.semi.systems`. -function DensityReinitializationCallback(system_indices::AbstractRange{<:Integer}; - interval::Integer=0, dt=0.0, +# Keywords +- `interval=0`: Reinitialize the density every `interval` time steps. +- `dt`: Reinitialize the density in regular intervals of `dt` in terms + of integration time. +- `reinit_initial_solution`: Reinitialize the initial solution. +""" +function DensityReinitializationCallback(system; interval::Integer=0, dt=0.0, reinit_initial_solution=true) - return reinit_callback_from_indices(system_indices; interval, dt, reinit_initial_solution) -end - -function reinit_callback_from_indices(system_indices; interval::Integer=0, dt=0.0, - reinit_initial_solution=true) if dt > 0 && interval > 0 error("Setting both interval and dt is not supported!") end @@ -99,43 +52,27 @@ function reinit_callback_from_indices(system_indices; interval::Integer=0, dt=0. interval = Float64(dt) end + check_density_reinit_system(system) + last_t = -Inf - system_indices_ = normalize_reinit_system_indices(system_indices) - reinit_cb = DensityReinitializationCallback(interval, nothing, system_indices_, - last_t, reinit_initial_solution) + reinit_cb = DensityReinitializationCallbackAffect(interval, system, last_t, + reinit_initial_solution) return DiscreteCallback(reinit_cb, reinit_cb, save_positions=(false, false), initialize=(initialize_reinit_cb!)) end -function DensityReinitializationCallback(; system_index=nothing, system_indices=nothing, - interval::Integer=0, - dt=0.0, reinit_initial_solution=true) - if !isnothing(system_index) && !isnothing(system_indices) - throw(ArgumentError("pass either `system_index` or `system_indices`, not both")) - end - - if isnothing(system_index) && isnothing(system_indices) - throw(ArgumentError("pass `system_index` or `system_indices`")) - end - - indices = isnothing(system_index) ? system_indices : system_index - - return reinit_callback_from_indices(indices; interval, dt, reinit_initial_solution) -end - function initialize_reinit_cb!(cb, u, t, integrator) initialize_reinit_cb!(cb.affect!, u, t, integrator) end -function initialize_reinit_cb!(cb::DensityReinitializationCallback, u, t, integrator) +function initialize_reinit_cb!(cb::DensityReinitializationCallbackAffect, u, t, integrator) semi = integrator.p.semi foreach_reinit_system(cb, semi) do system check_density_reinit_system(system) end - # Reinitialize initial solution if cb.reinit_initial_solution # Update systems to compute quantities like density and pressure. v_ode, u_ode = u.x @@ -151,21 +88,22 @@ function initialize_reinit_cb!(cb::DensityReinitializationCallback, u, t, integr end # condition with interval -function (reinit_callback::DensityReinitializationCallback{Int})(u, t, integrator) +function (reinit_callback::DensityReinitializationCallbackAffect{<:Integer})(u, t, + integrator) (; interval) = reinit_callback return condition_integrator_interval(integrator, interval, save_final_solution=false) end # condition with dt -function (reinit_callback::DensityReinitializationCallback)(u, t, integrator) +function (reinit_callback::DensityReinitializationCallbackAffect)(u, t, integrator) (; interval, last_t) = reinit_callback return (t - last_t) >= interval end # affect! -function (reinit_callback::DensityReinitializationCallback)(integrator) +function (reinit_callback::DensityReinitializationCallbackAffect)(integrator) vu_ode = integrator.u semi = integrator.p.semi @@ -180,7 +118,7 @@ function (reinit_callback::DensityReinitializationCallback)(integrator) return integrator end -function reinitialize_density!(reinit_callback::DensityReinitializationCallback, +function reinitialize_density!(reinit_callback::DensityReinitializationCallbackAffect, vu_ode, semi) v_ode, u_ode = vu_ode.x @@ -195,24 +133,15 @@ function reinitialize_density!(reinit_callback::DensityReinitializationCallback, return reinit_callback end -function foreach_reinit_system(f, reinit_callback::DensityReinitializationCallback, semi) - (; particle_system, system_indices) = reinit_callback - - if !isnothing(system_indices) - foreach(system_indices) do system_index - f(current_reinit_system(system_index, semi)) - end - - return nothing - end - - system_index = findfirst(system -> system === particle_system, semi.systems) +function foreach_reinit_system(f, reinit_callback::DensityReinitializationCallbackAffect, + semi) + system_index = findfirst(s -> s === reinit_callback.system, semi.systems) if isnothing(system_index) - throw(ArgumentError("system is not in the semidiscretization. This can happen " * - "when `semidiscretize` replaced the system. Pass a system " * - "index to `DensityReinitializationCallback` to select the " * - "system independently of replacement.")) + throw(ArgumentError("system is not in the semidiscretization. This can " * + "happen when `semidiscretize` creates a copy of the " * + "system. Create the callback with the corresponding " * + "system from `ode.p.semi.systems`.")) end f(semi.systems[system_index]) @@ -220,15 +149,6 @@ function foreach_reinit_system(f, reinit_callback::DensityReinitializationCallba return nothing end -function current_reinit_system(system_index, semi) - if system_index > length(semi.systems) - throw(ArgumentError("system index $system_index is out of bounds for a " * - "semidiscretization with $(length(semi.systems)) systems")) - end - - return semi.systems[system_index] -end - function check_density_reinit_system(particle_system) if !hasproperty(particle_system, :density_calculator) throw(ArgumentError("density reinitialization requires a system with a density calculator")) @@ -240,30 +160,3 @@ function check_density_reinit_system(particle_system) return particle_system end - -function normalize_reinit_system_indices(system_index::Integer) - return (normalize_reinit_system_index(system_index),) -end - -function normalize_reinit_system_indices(system_indices) - system_indices_ = Tuple(normalize_reinit_system_index(system_index) - for system_index in system_indices) - - if isempty(system_indices_) - throw(ArgumentError("`system_indices` must not be empty")) - end - - if !allunique(system_indices_) - throw(ArgumentError("`system_indices` must be unique")) - end - - return system_indices_ -end - -function normalize_reinit_system_index(system_index::Integer) - if system_index < 1 - throw(ArgumentError("system indices must be positive")) - end - - return Int(system_index) -end diff --git a/test/callbacks/density_reinit.jl b/test/callbacks/density_reinit.jl index 1587378569..3c9aff3163 100644 --- a/test/callbacks/density_reinit.jl +++ b/test/callbacks/density_reinit.jl @@ -8,6 +8,12 @@ name::Symbol end + struct MockDensityReinitIntegrator + p::Any + u::Any + t::Float64 + end + density_reinit_calls = Symbol[] TrixiParticles.wrap_v(v_ode, system::MockDensityReinitSystem, semi) = (:v, system.name) @@ -45,6 +51,43 @@ @test callback(nothing, 0.1, nothing) end + @testset verbose=true "reinit initial solution" begin + system = MockDensityReinitSystem(nothing, :fluid) + callback = DensityReinitializationCallback(system; interval=1, + reinit_initial_solution=false).affect! + + @test !callback.reinit_initial_solution + end + + @testset verbose=true "initialize respects reinit_initial_solution" begin + empty!(density_reinit_calls) + + system = MockDensityReinitSystem(nothing, :fluid) + vu_ode = (; x=(:v_ode, :u_ode)) + semi = (; systems=(system,)) + integrator = MockDensityReinitIntegrator((; semi), vu_ode, 0.0) + + TrixiParticles.get_neighborhood_search(system::MockDensityReinitSystem, + neighbor::MockDensityReinitSystem, + semi) = nothing + TrixiParticles.update_nhs!(neighborhood_search::Nothing, + system::MockDensityReinitSystem, + neighbor::MockDensityReinitSystem, + u_system, u_neighbor, semi) = nothing + TrixiParticles.u_modified!(integrator::MockDensityReinitIntegrator, + is_modified) = nothing + + callback = DensityReinitializationCallback(system; interval=1, + reinit_initial_solution=false).affect! + TrixiParticles.initialize_reinit_cb!(callback, vu_ode, 0.0, integrator) + @test isempty(density_reinit_calls) + + callback = DensityReinitializationCallback(system; interval=1, + reinit_initial_solution=true).affect! + TrixiParticles.initialize_reinit_cb!(callback, vu_ode, 0.0, integrator) + @test density_reinit_calls == [:fluid] + end + @testset verbose=true "selected system" begin empty!(density_reinit_calls) @@ -59,61 +102,40 @@ @test density_reinit_calls == [:fluid1] end - @testset verbose=true "selected system after semidiscretization replacement" begin + @testset verbose=true "selected semidiscretized system" begin empty!(density_reinit_calls) - replacement_system1 = MockDensityReinitSystem(Ref(:replacement1), :replacement1) - replacement_system2 = MockDensityReinitSystem(Ref(:replacement2), :replacement2) - - callback = DensityReinitializationCallback(system_index=1, interval=1).affect! + # Simulate the case where `semidiscretize` creates a copy of the system. + original_system = MockDensityReinitSystem(nothing, :original) + semidiscretized_system = MockDensityReinitSystem(nothing, :semidiscretized) vu_ode = (; x=(:v_ode, :u_ode)) - semi_replaced = (; systems=(replacement_system1, replacement_system2)) + semi = (; systems=(semidiscretized_system,)) - TrixiParticles.reinitialize_density!(callback, vu_ode, semi_replaced) + callback = DensityReinitializationCallback(semidiscretized_system; + interval=1).affect! + TrixiParticles.reinitialize_density!(callback, vu_ode, semi) + @test density_reinit_calls == [:semidiscretized] - @test density_reinit_calls == [:replacement1] + callback = DensityReinitializationCallback(original_system; interval=1).affect! + @test_throws ArgumentError TrixiParticles.reinitialize_density!(callback, vu_ode, + semi) end - @testset verbose=true "selected systems after semidiscretization replacement" begin + @testset verbose=true "system validation" begin empty!(density_reinit_calls) - replacement_system1 = MockDensityReinitSystem(Ref(:replacement1), :replacement1) - replacement_system2 = MockDensityReinitSystem(Ref(:replacement2), :replacement2) - replacement_system3 = MockDensityReinitSystem(Ref(:replacement3), :replacement3) - - callback = DensityReinitializationCallback(system_indices=(3, 1), - interval=1).affect! - vu_ode = (; x=(:v_ode, :u_ode)) - semi_replaced = (; systems=(replacement_system1, replacement_system2, - replacement_system3)) - - TrixiParticles.reinitialize_density!(callback, vu_ode, semi_replaced) - - @test density_reinit_calls == [:replacement3, :replacement1] - end - - @testset verbose=true "system index validation" begin - @test_throws ArgumentError DensityReinitializationCallback(system_index=0, - interval=1) - @test_throws ArgumentError DensityReinitializationCallback(system_indices=(), + system = MockDensityReinitSystem(nothing, :fluid) + @test_throws MethodError DensityReinitializationCallback(; interval=1) + @test_throws MethodError DensityReinitializationCallback(system, system; interval=1) + @test_throws ArgumentError DensityReinitializationCallback(MockDensityReinitSystem(SummationDensity(), + :fluid); interval=1) - @test_throws ArgumentError DensityReinitializationCallback(system_indices=(1, 1), + @test_throws ArgumentError DensityReinitializationCallback(MockNoDensityReinitSystem(:boundary); interval=1) - callback = DensityReinitializationCallback(system_index=3, interval=1).affect! + callback = DensityReinitializationCallback(system; interval=1).affect! vu_ode = (; x=(:v_ode, :u_ode)) - semi = (; systems=(MockDensityReinitSystem(nothing, :fluid),)) - - @test_throws ArgumentError TrixiParticles.reinitialize_density!(callback, vu_ode, - semi) - - callback = DensityReinitializationCallback(system_index=1, interval=1).affect! - semi = (; systems=(MockDensityReinitSystem(SummationDensity(), :fluid),)) - - @test_throws ArgumentError TrixiParticles.reinitialize_density!(callback, vu_ode, - semi) - - semi = (; systems=(MockNoDensityReinitSystem(:boundary),)) + semi = (; systems=(MockDensityReinitSystem(nothing, :other_fluid),)) @test_throws ArgumentError TrixiParticles.reinitialize_density!(callback, vu_ode, semi) From 0f0e9a27c1a2d2c89273756f4bbfec61a31a2528 Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Mon, 18 May 2026 13:22:38 +0200 Subject: [PATCH 5/5] simplify api --- examples/fluid/dam_break_2d.jl | 5 +-- src/callbacks/density_reinit.jl | 50 ++++++++++++--------------- test/callbacks/density_reinit.jl | 58 ++++++++++++++++++++------------ 3 files changed, 61 insertions(+), 52 deletions(-) diff --git a/examples/fluid/dam_break_2d.jl b/examples/fluid/dam_break_2d.jl index ee6f837dfa..11ee987599 100644 --- a/examples/fluid/dam_break_2d.jl +++ b/examples/fluid/dam_break_2d.jl @@ -111,9 +111,10 @@ extra_callback = nothing extra_callback2 = nothing use_reinit = false -# Use the system stored in `ode.p.semi`, since `semidiscretize` can replace systems. +# The callback stores the system index, so it remains valid if `semidiscretize` +# replaces systems internally. density_reinit_cb = use_reinit ? - DensityReinitializationCallback(ode.p.semi.systems[1]; interval=10) : + DensityReinitializationCallback(fluid_system, semi; interval=10) : nothing stepsize_callback = StepsizeCallback(cfl=0.9) diff --git a/src/callbacks/density_reinit.jl b/src/callbacks/density_reinit.jl index 788c263ca1..8de4f87fe3 100644 --- a/src/callbacks/density_reinit.jl +++ b/src/callbacks/density_reinit.jl @@ -1,6 +1,6 @@ -mutable struct DensityReinitializationCallbackAffect{I, PS} +mutable struct DensityReinitializationCallbackAffect{I} interval::I - system::PS + system_index::Int last_t::Float64 reinit_initial_solution::Bool end @@ -27,14 +27,14 @@ function Base.show(io::IO, ::MIME"text/plain", end """ - DensityReinitializationCallback(system; interval::Integer=0, dt=0.0, + DensityReinitializationCallback(system, semi; interval::Integer=0, dt=0.0, reinit_initial_solution=true) Callback to reinitialize the density field when using [`ContinuityDensity`](@ref) [Panizzo2007](@cite). -Pass a system to reinitialize its density field. The system must be the system stored -in the semidiscretization used by the integrator. If [`semidiscretize`](@ref) creates -a copy of the system, pass the corresponding system from `ode.p.semi.systems`. +Pass `system` and the [`Semidiscretization`](@ref) containing it. The callback stores +the system index and uses the corresponding system from the integrator semidiscretization +at runtime, which remains valid if [`semidiscretize`](@ref) replaces systems internally. # Keywords - `interval=0`: Reinitialize the density every `interval` time steps. @@ -44,7 +44,7 @@ a copy of the system, pass the corresponding system from `ode.p.semi.systems`. the first solver step after each `dt` interval has elapsed. - `reinit_initial_solution`: Reinitialize the initial solution (default=true). """ -function DensityReinitializationCallback(system; interval::Integer=0, dt=0.0, +function DensityReinitializationCallback(system, semi; interval::Integer=0, dt=0.0, reinit_initial_solution=true) if dt > 0 && interval > 0 error("Setting both interval and dt is not supported!") @@ -56,9 +56,10 @@ function DensityReinitializationCallback(system; interval::Integer=0, dt=0.0, check_density_reinit_system(system) + system_index = system_indices(system, semi) last_t = -Inf - reinit_cb = DensityReinitializationCallbackAffect(interval, system, last_t, + reinit_cb = DensityReinitializationCallbackAffect(interval, system_index, last_t, reinit_initial_solution) return DiscreteCallback(reinit_cb, reinit_cb, save_positions=(false, false), @@ -71,9 +72,7 @@ end function initialize_reinit_cb!(cb::DensityReinitializationCallbackAffect, u, t, integrator) semi = integrator.p.semi - foreach_reinit_system(cb, semi) do system - check_density_reinit_system(system) - end + check_density_reinit_system(current_reinit_system(cb.system_index, semi)) if cb.reinit_initial_solution # Update systems to compute quantities like density and pressure. @@ -124,31 +123,24 @@ function reinitialize_density!(reinit_callback::DensityReinitializationCallbackA vu_ode, semi) v_ode, u_ode = vu_ode.x - foreach_reinit_system(reinit_callback, semi) do particle_system - check_density_reinit_system(particle_system) - v = wrap_v(v_ode, particle_system, semi) - u = wrap_u(u_ode, particle_system, semi) + particle_system = current_reinit_system(reinit_callback.system_index, semi) + check_density_reinit_system(particle_system) - reinit_density!(particle_system, v, u, v_ode, u_ode, semi) - end + v = wrap_v(v_ode, particle_system, semi) + u = wrap_u(u_ode, particle_system, semi) + + reinit_density!(particle_system, v, u, v_ode, u_ode, semi) return reinit_callback end -function foreach_reinit_system(f, reinit_callback::DensityReinitializationCallbackAffect, - semi) - system_index = findfirst(s -> s === reinit_callback.system, semi.systems) - - if isnothing(system_index) - throw(ArgumentError("system is not in the semidiscretization. This can " * - "happen when `semidiscretize` creates a copy of the " * - "system. Create the callback with the corresponding " * - "system from `ode.p.semi.systems`.")) +function current_reinit_system(system_index, semi) + if !(1 <= system_index <= length(semi.systems)) + throw(ArgumentError("system index $system_index is out of bounds for a " * + "semidiscretization with $(length(semi.systems)) systems")) end - f(semi.systems[system_index]) - - return nothing + return semi.systems[system_index] end function check_density_reinit_system(particle_system) diff --git a/test/callbacks/density_reinit.jl b/test/callbacks/density_reinit.jl index 473745b3a2..a9d0386593 100644 --- a/test/callbacks/density_reinit.jl +++ b/test/callbacks/density_reinit.jl @@ -27,7 +27,8 @@ @testset verbose=true "show" begin system = MockDensityReinitSystem(nothing, :fluid) - callback = DensityReinitializationCallback(system; interval=10) + semi = (; systems=(system,)) + callback = DensityReinitializationCallback(system, semi; interval=10) show_compact = "DensityReinitializationCallback(interval=10, reinit_initial_solution=true)" @test repr(callback) == show_compact @@ -44,7 +45,8 @@ @testset verbose=true "dt condition" begin system = MockDensityReinitSystem(nothing, :fluid) - callback = DensityReinitializationCallback(system; dt=0.1).affect! + semi = (; systems=(system,)) + callback = DensityReinitializationCallback(system, semi; dt=0.1).affect! callback.last_t = 0.0 @test !callback(nothing, 0.099, nothing) @@ -54,7 +56,8 @@ @testset verbose=true "reinit initial solution" begin system = MockDensityReinitSystem(nothing, :fluid) - callback = DensityReinitializationCallback(system; interval=1, + semi = (; systems=(system,)) + callback = DensityReinitializationCallback(system, semi; interval=1, reinit_initial_solution=false).affect! @test !callback.reinit_initial_solution @@ -78,12 +81,12 @@ TrixiParticles.u_modified!(integrator::MockDensityReinitIntegrator, is_modified) = nothing - callback = DensityReinitializationCallback(system; interval=1, + callback = DensityReinitializationCallback(system, semi; interval=1, reinit_initial_solution=false).affect! TrixiParticles.initialize_reinit_cb!(callback, vu_ode, 0.0, integrator) @test isempty(density_reinit_calls) - callback = DensityReinitializationCallback(system; interval=1, + callback = DensityReinitializationCallback(system, semi; interval=1, reinit_initial_solution=true).affect! TrixiParticles.initialize_reinit_cb!(callback, vu_ode, 0.0, integrator) @test density_reinit_calls == [:fluid] @@ -94,9 +97,11 @@ system1 = MockDensityReinitSystem(nothing, :fluid1) system2 = MockDensityReinitSystem(nothing, :fluid2) - callback = DensityReinitializationCallback(system1; interval=1).affect! vu_ode = (; x=(:v_ode, :u_ode)) semi = (; systems=(system1, system2)) + callback = DensityReinitializationCallback(system1, semi; interval=1).affect! + + @test callback.system_index == 1 TrixiParticles.reinitialize_density!(callback, vu_ode, semi) @@ -106,39 +111,50 @@ @testset verbose=true "selected semidiscretized system" begin empty!(density_reinit_calls) - # Simulate the case where `semidiscretize` creates a copy of the system. - original_system = MockDensityReinitSystem(nothing, :original) - semidiscretized_system = MockDensityReinitSystem(nothing, :semidiscretized) + original_system1 = MockDensityReinitSystem(nothing, :original1) + original_system2 = MockDensityReinitSystem(nothing, :original2) + semidiscretized_system1 = MockDensityReinitSystem(nothing, :semidiscretized1) + semidiscretized_system2 = MockDensityReinitSystem(nothing, :semidiscretized2) vu_ode = (; x=(:v_ode, :u_ode)) - semi = (; systems=(semidiscretized_system,)) + semi = (; systems=(original_system1, original_system2)) + semi_replaced = (; systems=(semidiscretized_system1, semidiscretized_system2)) - callback = DensityReinitializationCallback(semidiscretized_system; + # Simulate the case where `semidiscretize` creates a copy of the system. + callback = DensityReinitializationCallback(original_system2, semi; interval=1).affect! - TrixiParticles.reinitialize_density!(callback, vu_ode, semi) - @test density_reinit_calls == [:semidiscretized] + @test callback.system_index == 2 + + TrixiParticles.reinitialize_density!(callback, vu_ode, semi_replaced) + @test density_reinit_calls == [:semidiscretized2] - callback = DensityReinitializationCallback(original_system; interval=1).affect! @test_throws ArgumentError TrixiParticles.reinitialize_density!(callback, vu_ode, - semi) + (; systems=())) end @testset verbose=true "system validation" begin empty!(density_reinit_calls) system = MockDensityReinitSystem(nothing, :fluid) + semi = (; systems=(system,)) + other_semi = (; systems=(MockDensityReinitSystem(nothing, :other_fluid),)) + @test_throws MethodError DensityReinitializationCallback(; interval=1) - @test_throws MethodError DensityReinitializationCallback(system, system; interval=1) + @test_throws MethodError DensityReinitializationCallback(system; interval=1) @test_throws ArgumentError DensityReinitializationCallback(MockDensityReinitSystem(SummationDensity(), - :fluid); + :fluid), + semi; + interval=1) + @test_throws ArgumentError DensityReinitializationCallback(MockNoDensityReinitSystem(:boundary), + semi; interval=1) - @test_throws ArgumentError DensityReinitializationCallback(MockNoDensityReinitSystem(:boundary); + @test_throws ArgumentError DensityReinitializationCallback(system, other_semi; interval=1) - callback = DensityReinitializationCallback(system; interval=1).affect! + callback = DensityReinitializationCallback(system, semi; interval=1).affect! vu_ode = (; x=(:v_ode, :u_ode)) - semi = (; systems=(MockDensityReinitSystem(nothing, :other_fluid),)) @test_throws ArgumentError TrixiParticles.reinitialize_density!(callback, vu_ode, - semi) + (; + systems=(MockNoDensityReinitSystem(:boundary),))) end end