From f8c99f7d85923a9ade5c8e7fcea57df908c47236 Mon Sep 17 00:00:00 2001
From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com>
Date: Mon, 21 Apr 2025 20:58:17 +0200
Subject: [PATCH 1/9] Implement `SplitIntegrationCallback` to integrate TLSPH
 separately

---
 src/TrixiParticles.jl                         |   7 +-
 src/callbacks/callbacks.jl                    |   1 +
 src/callbacks/split_integration.jl            | 242 ++++++++++++++++++
 src/general/semidiscretization.jl             | 192 +++++++++-----
 src/schemes/solid/total_lagrangian_sph/rhs.jl |  15 +-
 5 files changed, 384 insertions(+), 73 deletions(-)
 create mode 100644 src/callbacks/split_integration.jl

diff --git a/src/TrixiParticles.jl b/src/TrixiParticles.jl
index 94e4a639b2..98a6f87c34 100644
--- a/src/TrixiParticles.jl
+++ b/src/TrixiParticles.jl
@@ -21,8 +21,9 @@ using Polyester: Polyester, @batch
 using Printf: @printf, @sprintf
 using RecipesBase: RecipesBase, @series
 using Random: seed!
-using SciMLBase: CallbackSet, DiscreteCallback, DynamicalODEProblem, u_modified!,
-                 get_tmp_cache, set_proposed_dt!, ODESolution, ODEProblem, terminate!
+using SciMLBase: SciMLBase, CallbackSet, DiscreteCallback, DynamicalODEProblem, u_modified!,
+                 get_tmp_cache, set_proposed_dt!, ODESolution, ODEProblem, terminate!,
+                 add_tstop!
 @reexport using StaticArrays: SVector
 using StaticArrays: @SMatrix, SMatrix, setindex
 using StrideArrays: PtrArray, StaticInt
@@ -62,7 +63,7 @@ export WeaklyCompressibleSPHSystem, EntropicallyDampedSPHSystem, TotalLagrangian
 export BoundaryZone, InFlow, OutFlow, BidirectionalFlow
 export InfoCallback, SolutionSavingCallback, DensityReinitializationCallback,
        PostprocessCallback, StepsizeCallback, UpdateCallback, SteadyStateReachedCallback,
-       ParticleShiftingCallback
+       ParticleShiftingCallback, SplitIntegrationCallback
 export ContinuityDensity, SummationDensity
 export PenaltyForceGanzenmueller, TransportVelocityAdami
 export SchoenbergCubicSplineKernel, SchoenbergQuarticSplineKernel,
diff --git a/src/callbacks/callbacks.jl b/src/callbacks/callbacks.jl
index 062a16e3a9..5d058794fa 100644
--- a/src/callbacks/callbacks.jl
+++ b/src/callbacks/callbacks.jl
@@ -33,3 +33,4 @@ include("stepsize.jl")
 include("update.jl")
 include("steady_state_reached.jl")
 include("particle_shifting.jl")
+include("split_integration.jl")
diff --git a/src/callbacks/split_integration.jl b/src/callbacks/split_integration.jl
new file mode 100644
index 0000000000..adee792758
--- /dev/null
+++ b/src/callbacks/split_integration.jl
@@ -0,0 +1,242 @@
+mutable struct SplitIntegrationCallback
+    integrator :: Any
+    alg        :: Any
+    kwargs     :: Any
+end
+
+@doc raw"""
+    SplitIntegrationCallback()
+
+Callback to...
+"""
+function SplitIntegrationCallback(alg; kwargs...)
+    split_integration_callback = SplitIntegrationCallback(nothing, alg, kwargs)
+
+    # The first one is the `condition`, the second the `affect!`
+    return DiscreteCallback(split_integration_callback, split_integration_callback,
+                            initialize=initialize_split_integration!,
+                            save_positions=(false, false))
+end
+
+function initialize_split_integration!(cb, u, t, integrator)
+    semi = integrator.p
+    split_integration_callback = cb.affect!
+    (; alg, kwargs) = split_integration_callback
+
+    # Disable TLSPH integration in the original integrator
+    semi.integrate_tlsph[] = false
+
+    # Create split integrator with TLSPH systems only
+    systems = filter(i -> i isa TotalLagrangianSPHSystem, semi.systems)
+
+    # These neighborhood searches are never used
+    semi_split = Semidiscretization(systems...,
+                                    neighborhood_search=TrivialNeighborhoodSearch{ndims(first(systems))}(),
+                                    parallelization_backend=semi.parallelization_backend)
+
+    sizes_u = (u_nvariables(system) * n_moving_particles(system) for system in systems)
+    sizes_v = (v_nvariables(system) * n_moving_particles(system) for system in systems)
+
+    v_ode, u_ode = integrator.u.x
+    v0_ode_split = similar(v_ode, sum(sizes_v))
+    u0_ode_split = similar(u_ode, sum(sizes_u))
+
+    # Copy the initial conditions to the split integrator
+    foreach_system(semi_split) do system
+        v = wrap_v(v_ode, system, semi)
+        u = wrap_u(u_ode, system, semi)
+        v_split = wrap_v(v0_ode_split, system, semi_split)
+        u_split = wrap_u(u0_ode_split, system, semi_split)
+
+        @threaded semi for particle in each_moving_particle(system)
+            for i in axes(v, 1)
+                v_split[i, particle] = v[i, particle]
+            end
+
+            for i in axes(u, 1)
+                u_split[i, particle] = u[i, particle]
+            end
+        end
+    end
+
+    # No tstops by default
+    tspan = (integrator.t, integrator.t + 1)
+    p = (; v_ode, u_ode, semi, semi_split)
+    ode_split = DynamicalODEProblem(kick_split!, drift_split!, v0_ode_split, u0_ode_split,
+                                    tspan, p)
+
+    # Create the split integrator
+    split_integrator = SciMLBase.init(ode_split, alg; save_everystep=false, kwargs...)
+    # A zero `tspan` sets `tdir` to zero, which breaks adding tstops
+    SciMLBase.pop_tstop!(split_integrator)
+    split_integration_callback.integrator = split_integrator
+
+    return cb
+end
+
+# `condition`
+function (split_integration_callback::SplitIntegrationCallback)(u, t, integrator)
+    return true
+end
+
+# `affect!`
+function (split_integration_callback::SplitIntegrationCallback)(integrator)
+    split_integrator = split_integration_callback.integrator
+    semi_split = split_integrator.p.semi_split
+
+    semi = integrator.p
+    new_t = integrator.t
+
+    v_ode, u_ode = integrator.u.x
+
+    @assert semi == split_integrator.p.semi
+    split_integrator.p = (; v_ode, u_ode, semi, semi_split)
+
+    @trixi_timeit timer() "split integration" 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.
+        update_systems_and_nhs(v_ode, u_ode, semi, new_t; update_from_callback=true)
+
+        # Integrate the split integrator to the new time
+        add_tstop!(split_integrator, new_t)
+        SciMLBase.solve!(split_integrator)
+
+        v_ode_split, u_ode_split = split_integrator.u.x
+
+        # Copy the solutions back to the original integrator
+        foreach_system(semi_split) do system
+            v = wrap_v(v_ode, system, semi)
+            u = wrap_u(u_ode, system, semi)
+            v_split = wrap_v(v_ode_split, system, semi_split)
+            u_split = wrap_u(u_ode_split, system, semi_split)
+
+            @threaded semi for particle in each_moving_particle(system)
+                for i in axes(v, 1)
+                    v[i, particle] = v_split[i, particle]
+                end
+
+                for i in axes(u, 1)
+                    u[i, particle] = u_split[i, particle]
+                end
+            end
+        end
+    end
+
+    # Tell OrdinaryDiffEq that `u` has been modified
+    u_modified!(integrator, true)
+
+    return integrator
+end
+
+function kick_split!(dv_ode_split, v_ode_split, u_ode_split, p, t)
+    (; v_ode, u_ode, semi, semi_split) = p
+
+    @trixi_timeit timer() "split integration kick!" begin
+        @trixi_timeit timer() "reset ∂v/∂t" set_zero!(dv_ode_split)
+
+        # Only loop over TLSPH systems
+        @trixi_timeit timer() "copy to large v,u" foreach_system(semi_split) do system
+            v = wrap_v(v_ode, system, semi)
+            u = wrap_u(u_ode, system, semi)
+            v_split = wrap_v(v_ode_split, system, semi_split)
+            u_split = wrap_u(u_ode_split, system, semi_split)
+
+            @threaded semi for particle in each_moving_particle(system)
+                for i in axes(v, 1)
+                    v[i, particle] = v_split[i, particle]
+                end
+
+                for i in axes(u, 1)
+                    u[i, particle] = u_split[i, particle]
+                end
+            end
+        end
+
+        update_nhs_fun = (semi, u_ode) -> update_nhs_split!(semi, u_ode, u_ode_split,
+                                                            semi_split)
+        @trixi_timeit timer() "update systems and nhs" update_systems_and_nhs(v_ode, u_ode,
+                                                                              semi, t;
+                                                                              systems=semi_split.systems,
+                                                                              update_nhs_fun)
+
+        @trixi_timeit timer() "system interaction" begin
+            # Only loop over systems in the split integrator
+            foreach_system(semi_split) do system
+                # Loop over neighbors in the big integrator
+                foreach_system(semi) do neighbor
+                    # Construct string for the interactions timer.
+                    # Avoid allocations from string construction when no timers are used.
+                    if timeit_debug_enabled()
+                        system_index = system_indices(system, semi)
+                        neighbor_index = system_indices(neighbor, semi)
+                        timer_str = "$(timer_name(system))$system_index-$(timer_name(neighbor))$neighbor_index"
+                    else
+                        timer_str = ""
+                    end
+
+                    dv = wrap_v(dv_ode_split, system, semi_split)
+                    v_system = wrap_v(v_ode_split, system, semi_split)
+                    u_system = wrap_u(u_ode_split, system, semi_split)
+
+                    v_neighbor = wrap_v(v_ode, neighbor, semi)
+                    u_neighbor = wrap_u(u_ode, neighbor, semi)
+
+                    @trixi_timeit timer() timer_str begin
+                        interact!(dv, v_system, u_system, v_neighbor, u_neighbor,
+                                  system, neighbor, semi; integrate_tlsph=true)
+                    end
+                end
+            end
+        end
+
+        @trixi_timeit timer() "source terms" add_source_terms!(dv_ode_split, v_ode_split,
+                                                               u_ode_split, semi, t;
+                                                               semi_wrap=semi_split)
+    end
+end
+
+function drift_split!(du_ode, v_ode, u_ode, p, t)
+    drift!(du_ode, v_ode, u_ode, p.semi_split, t)
+end
+
+function update_nhs_split!(semi, u_ode, u_ode_split, semi_split)
+    # Only loop over systems in the split integrator
+    foreach_system(semi_split) do system
+        u_system = wrap_u(u_ode_split, system, semi_split)
+
+        # Loop over neighbors in the big integrator
+        foreach_system(semi) do neighbor
+            # Static NHS for solid-solid (same system) and no interaction for two distinct
+            # solid systems (TODO).
+            if !(neighbor isa TotalLagrangianSPHSystem)
+                u_neighbor = wrap_u(u_ode, neighbor, semi)
+                neighborhood_search = get_neighborhood_search(system, neighbor, semi)
+
+                # Only the TLSPH particles are moving. All other systems are frozen.
+                # Note that this does nothing when a grid NHS is used.
+                update_nhs!(neighborhood_search, system, neighbor, u_system, u_neighbor,
+                            semi; points_moving=(true, false))
+            end
+        end
+    end
+end
+
+function Base.show(io::IO, cb::DiscreteCallback{<:Any, <:SplitIntegrationCallback})
+    @nospecialize cb # reduce precompilation time
+    print(io, "SplitIntegrationCallback(integrator=", cb.affect!.integrator, ")")
+end
+
+function Base.show(io::IO, ::MIME"text/plain",
+                   cb::DiscreteCallback{<:Any, <:SplitIntegrationCallback})
+    @nospecialize cb # reduce precompilation time
+
+    if get(io, :compact, false)
+        show(io, cb)
+    else
+        update_cb = cb.affect!
+        setup = [
+            "integrator" => update_cb.integrator
+        ]
+        summary_box(io, "SplitIntegrationCallback", setup)
+    end
+end
diff --git a/src/general/semidiscretization.jl b/src/general/semidiscretization.jl
index c0a70da3f3..57724d9763 100644
--- a/src/general/semidiscretization.jl
+++ b/src/general/semidiscretization.jl
@@ -47,23 +47,24 @@ semi = Semidiscretization(fluid_system, boundary_system,
 └──────────────────────────────────────────────────────────────────────────────────────────────────┘
 ```
 """
-struct Semidiscretization{BACKEND, S, RU, RV, NS}
+struct Semidiscretization{BACKEND, S, RU, RV, NS, IT}
     systems                 :: S
     ranges_u                :: RU
     ranges_v                :: RV
     neighborhood_searches   :: NS
     parallelization_backend :: BACKEND
+    integrate_tlsph         :: IT
 
     # Dispatch at `systems` to distinguish this constructor from the one below when
     # 4 systems are passed.
     # This is an internal constructor only used in `test/count_allocations.jl`
     # and by Adapt.jl.
     function Semidiscretization(systems::Tuple, ranges_u, ranges_v, neighborhood_searches,
-                                parallelization_backend)
+                                parallelization_backend, integrate_tlsph)
         new{typeof(parallelization_backend), typeof(systems), typeof(ranges_u),
-            typeof(ranges_v), typeof(neighborhood_searches)}(systems, ranges_u, ranges_v,
-                                                             neighborhood_searches,
-                                                             parallelization_backend)
+            typeof(ranges_v), typeof(neighborhood_searches),
+            typeof(integrate_tlsph)}(systems, ranges_u, ranges_v, neighborhood_searches,
+                                     parallelization_backend, integrate_tlsph)
     end
 end
 
@@ -93,7 +94,7 @@ function Semidiscretization(systems...;
                      for system in systems)
 
     return Semidiscretization(systems, ranges_u, ranges_v, searches,
-                              parallelization_backend)
+                              parallelization_backend, Ref(true))
 end
 
 # Inline show function e.g. Semidiscretization(neighborhood_search=...)
@@ -431,7 +432,7 @@ function drift!(du_ode, v_ode, u_ode, semi, t)
 
                 @threaded semi for particle in each_moving_particle(system)
                     # This can be dispatched per system
-                    add_velocity!(du, v, particle, system)
+                    add_velocity!(du, v, particle, system, semi)
                 end
             end
         end
@@ -440,6 +441,18 @@ function drift!(du_ode, v_ode, u_ode, semi, t)
     return du_ode
 end
 
+@inline function add_velocity!(du, v, particle, system, semi::Semidiscretization)
+    add_velocity!(du, v, particle, system)
+end
+
+@inline function add_velocity!(du, v, particle, system::TotalLagrangianSPHSystem,
+                               semi::Semidiscretization)
+    # Only add velocity for TLSPH systems if they are integrated
+    if semi.integrate_tlsph[]
+        add_velocity!(du, v, particle, system)
+    end
+end
+
 @inline function add_velocity!(du, v, particle, system)
     for i in 1:ndims(system)
         du[i, particle] = v[i, particle]
@@ -467,11 +480,15 @@ function kick!(dv_ode, v_ode, u_ode, semi, t)
     return dv_ode
 end
 
-# Update the systems and neighborhood searches (NHS) for a simulation before calling `interact!` to compute forces
-function update_systems_and_nhs(v_ode, u_ode, semi, t; update_from_callback=false)
+# Update the systems and neighborhood searches (NHS) for a simulation
+# before calling `interact!` to compute forces.
+# `semi` and `update_nhs_fun` are overwritten by the `SplitIntegrationCallback`.
+function update_systems_and_nhs(v_ode, u_ode, semi, t;
+                                systems=semi, update_nhs_fun=update_nhs!,
+                                update_from_callback=false)
     # First update step before updating the NHS
     # (for example for writing the current coordinates in the solid system)
-    foreach_system(semi) do system
+    foreach_system(systems) do system
         v = wrap_v(v_ode, system, semi)
         u = wrap_u(u_ode, system, semi)
 
@@ -479,13 +496,13 @@ function update_systems_and_nhs(v_ode, u_ode, semi, t; update_from_callback=fals
     end
 
     # Update NHS
-    @trixi_timeit timer() "update nhs" update_nhs!(semi, u_ode)
+    @trixi_timeit timer() "update nhs" update_nhs_fun(semi, u_ode)
 
     # Second update step.
     # This is used to calculate density and pressure of the fluid systems
     # before updating the boundary systems,
     # since the fluid pressure is needed by the Adami interpolation.
-    foreach_system(semi) do system
+    foreach_system(systems) do system
         v = wrap_v(v_ode, system, semi)
         u = wrap_u(u_ode, system, semi)
 
@@ -493,7 +510,7 @@ function update_systems_and_nhs(v_ode, u_ode, semi, t; update_from_callback=fals
     end
 
     # Perform correction and pressure calculation
-    foreach_system(semi) do system
+    foreach_system(systems) do system
         v = wrap_v(v_ode, system, semi)
         u = wrap_u(u_ode, system, semi)
 
@@ -501,7 +518,7 @@ function update_systems_and_nhs(v_ode, u_ode, semi, t; update_from_callback=fals
     end
 
     # Final update step for all remaining systems
-    foreach_system(semi) do system
+    foreach_system(systems) do system
         v = wrap_v(v_ode, system, semi)
         u = wrap_u(u_ode, system, semi)
 
@@ -523,16 +540,22 @@ function update_nhs!(semi, u_ode)
     end
 end
 
-function add_source_terms!(dv_ode, v_ode, u_ode, semi, t)
-    foreach_system(semi) do system
-        dv = wrap_v(dv_ode, system, semi)
-        v = wrap_v(v_ode, system, semi)
-        u = wrap_u(u_ode, system, semi)
+# The `SplitIntegrationCallback` overwrites `semi_wrap` to use a different semi for wrapping
+# TODO `semi` is not used yet, but will be used when the source terms API is modified
+# to match the custom quantities API.
+function add_source_terms!(dv_ode, v_ode, u_ode, semi, t; semi_wrap=semi)
+    foreach_system(semi_wrap) do system
+        dv = wrap_v(dv_ode, system, semi_wrap)
+        v = wrap_v(v_ode, system, semi_wrap)
+        u = wrap_u(u_ode, system, semi_wrap)
 
         @threaded semi for particle in each_moving_particle(system)
-            # Dispatch by system type to exclude boundary systems
-            add_acceleration!(dv, particle, system)
-            add_source_terms_inner!(dv, v, u, particle, system, source_terms(system), t)
+            # Dispatch by system type to exclude boundary systems.
+            # `integrate_tlsph` is extracted from the `semi_wrap`, so that this function
+            # can be used in the `SplitIntegrationCallback` as well.
+            add_acceleration!(dv, particle, system, semi_wrap.integrate_tlsph[])
+            add_source_terms_inner!(dv, v, u, particle, system, source_terms(system), t,
+                                    semi_wrap.integrate_tlsph[])
         end
     end
 
@@ -542,6 +565,15 @@ end
 @inline source_terms(system) = nothing
 @inline source_terms(system::Union{FluidSystem, SolidSystem}) = system.source_terms
 
+@inline function add_acceleration!(dv, particle, system, integrate_tlsph)
+    add_acceleration!(dv, particle, system)
+end
+
+@inline function add_acceleration!(dv, particle, system::TotalLagrangianSPHSystem,
+                                   integrate_tlsph)
+    integrate_tlsph && add_acceleration!(dv, particle, system)
+end
+
 @inline add_acceleration!(dv, particle, system) = dv
 
 @inline function add_acceleration!(dv, particle, system::Union{FluidSystem, SolidSystem})
@@ -554,6 +586,17 @@ end
     return dv
 end
 
+@inline function add_source_terms_inner!(dv, v, u, particle, system, source_terms_, t,
+                                         integrate_tlsph)
+    add_source_terms_inner!(dv, v, u, particle, system, source_terms_, t)
+end
+
+@inline function add_source_terms_inner!(dv, v, u, particle,
+                                         system::TotalLagrangianSPHSystem,
+                                         source_terms_, t, integrate_tlsph)
+    integrate_tlsph && add_source_terms_inner!(dv, v, u, particle, system, source_terms_, t)
+end
+
 @inline function add_source_terms_inner!(dv, v, u, particle, system, source_terms_, t)
     coords = current_coords(u, system, particle)
     velocity = current_velocity(v, system, particle)
@@ -651,162 +694,175 @@ end
 function update_nhs!(neighborhood_search,
                      system::FluidSystem,
                      neighbor::Union{FluidSystem, TotalLagrangianSPHSystem},
-                     u_system, u_neighbor, semi)
-    # The current coordinates of fluids and solids change over time
+                     u_system, u_neighbor, semi;
+                     # The current coordinates of fluids and solids change over time
+                     points_moving=(true, true))
     update!(neighborhood_search,
             current_coordinates(u_system, system),
             current_coordinates(u_neighbor, neighbor),
-            semi, points_moving=(true, true))
+            semi; points_moving)
 end
 
 function update_nhs!(neighborhood_search,
                      system::FluidSystem, neighbor::BoundarySPHSystem,
-                     u_system, u_neighbor, semi)
-    # Boundary coordinates only change over time when `neighbor.ismoving[]`
+                     u_system, u_neighbor, semi;
+                     # Boundary coordinates only change over time when `neighbor.ismoving[]`
+                     points_moving=(true, neighbor.ismoving[]))
     update!(neighborhood_search,
             current_coordinates(u_system, system),
             current_coordinates(u_neighbor, neighbor),
-            semi, points_moving=(true, neighbor.ismoving[]))
+            semi; points_moving)
 end
 
 function update_nhs!(neighborhood_search,
                      system::FluidSystem, neighbor::OpenBoundarySPHSystem,
-                     u_system, u_neighbor, semi)
-    # The current coordinates of fluids and open boundaries change over time.
-
+                     u_system, u_neighbor, semi;
+                     # The current coordinates of fluids and open boundaries change over time
+                     points_moving=(true, true))
     # TODO: Update only `active_coordinates` of open boundaries.
     # Problem: Removing inactive particles from neighboring lists is necessary.
     update!(neighborhood_search,
             current_coordinates(u_system, system),
             current_coordinates(u_neighbor, neighbor),
-            semi, points_moving=(true, true))
+            semi; points_moving)
 end
 
 function update_nhs!(neighborhood_search,
                      system::OpenBoundarySPHSystem, neighbor::FluidSystem,
-                     u_system, u_neighbor, semi)
-    # The current coordinates of both open boundaries and fluids change over time.
-
+                     u_system, u_neighbor, semi;
+                     # The current coordinates of both open boundaries and fluids change over time
+                     points_moving=(true, true))
     # TODO: Update only `active_coordinates` of open boundaries.
     # Problem: Removing inactive particles from neighboring lists is necessary.
     update!(neighborhood_search,
             current_coordinates(u_system, system),
             current_coordinates(u_neighbor, neighbor),
-            semi, points_moving=(true, true))
+            semi; points_moving)
 end
 
 function update_nhs!(neighborhood_search,
                      system::OpenBoundarySPHSystem, neighbor::TotalLagrangianSPHSystem,
-                     u_system, u_neighbor)
+                     u_system, u_neighbor;
+                     points_moving=(false, false))
     # Don't update. This NHS is never used.
     return neighborhood_search
 end
 
 function update_nhs!(neighborhood_search,
                      system::TotalLagrangianSPHSystem, neighbor::OpenBoundarySPHSystem,
-                     u_system, u_neighbor)
+                     u_system, u_neighbor;
+                     points_moving=(false, false))
     # Don't update. This NHS is never used.
     return neighborhood_search
 end
 
 function update_nhs!(neighborhood_search,
                      system::TotalLagrangianSPHSystem, neighbor::FluidSystem,
-                     u_system, u_neighbor, semi)
+                     u_system, u_neighbor, semi;
+                     # The current coordinates of fluids and solids change over time
+                     points_moving=(true, true))
     # The current coordinates of fluids and solids change over time
     update!(neighborhood_search,
             current_coordinates(u_system, system),
             current_coordinates(u_neighbor, neighbor),
-            semi, points_moving=(true, true))
+            semi; points_moving)
 end
 
 function update_nhs!(neighborhood_search,
                      system::TotalLagrangianSPHSystem, neighbor::TotalLagrangianSPHSystem,
-                     u_system, u_neighbor, semi)
+                     u_system, u_neighbor, semi;
+                     points_moving=(false, false))
     # Don't update. Neighborhood search works on the initial coordinates, which don't change.
     return neighborhood_search
 end
 
 function update_nhs!(neighborhood_search,
                      system::TotalLagrangianSPHSystem, neighbor::BoundarySPHSystem,
-                     u_system, u_neighbor, semi)
+                     u_system, u_neighbor, semi;
+                     # The current coordinates of solids change over time
+                     points_moving=(true, neighbor.ismoving[]))
     # The current coordinates of solids change over time.
     # Boundary coordinates only change over time when `neighbor.ismoving[]`.
     update!(neighborhood_search,
             current_coordinates(u_system, system),
             current_coordinates(u_neighbor, neighbor),
-            semi, points_moving=(true, neighbor.ismoving[]))
+            semi; points_moving)
 end
 
 # This function is the same as the one below to avoid ambiguous dispatch when using `Union`
 function update_nhs!(neighborhood_search,
                      system::BoundarySPHSystem{<:BoundaryModelDummyParticles},
-                     neighbor::FluidSystem, u_system, u_neighbor, semi)
+                     neighbor::FluidSystem, u_system, u_neighbor, semi;
+                     # Boundary coordinates only change over time when `neighbor.ismoving[]`.
+                     # The current coordinates of fluids and solids change over time.
+                     points_moving=(system.ismoving[], true))
     # Depending on the density calculator of the boundary model, this NHS is used for
     # - kernel summation (`SummationDensity`)
     # - continuity equation (`ContinuityDensity`)
     # - pressure extrapolation (`AdamiPressureExtrapolation`)
-    #
-    # Boundary coordinates only change over time when `neighbor.ismoving[]`.
-    # The current coordinates of fluids and solids change over time.
     update!(neighborhood_search,
             current_coordinates(u_system, system),
             current_coordinates(u_neighbor, neighbor),
-            semi, points_moving=(system.ismoving[], true))
+            semi; points_moving)
 end
 
 # This function is the same as the one above to avoid ambiguous dispatch when using `Union`
 function update_nhs!(neighborhood_search,
                      system::BoundarySPHSystem{<:BoundaryModelDummyParticles},
-                     neighbor::TotalLagrangianSPHSystem, u_system, u_neighbor, semi)
+                     neighbor::TotalLagrangianSPHSystem, u_system, u_neighbor, semi;
+                     # Boundary coordinates only change over time when `neighbor.ismoving[]`.
+                     # The current coordinates of solids change over time.
+                     points_moving=(system.ismoving[], true))
     # Depending on the density calculator of the boundary model, this NHS is used for
     # - kernel summation (`SummationDensity`)
     # - continuity equation (`ContinuityDensity`)
     # - pressure extrapolation (`AdamiPressureExtrapolation`)
-    #
-    # Boundary coordinates only change over time when `neighbor.ismoving[]`.
-    # The current coordinates of fluids and solids change over time.
     update!(neighborhood_search,
             current_coordinates(u_system, system),
             current_coordinates(u_neighbor, neighbor),
-            semi, points_moving=(system.ismoving[], true))
+            semi; points_moving)
 end
 
 function update_nhs!(neighborhood_search,
                      system::BoundarySPHSystem{<:BoundaryModelDummyParticles},
                      neighbor::BoundarySPHSystem,
-                     u_system, u_neighbor, semi)
-    # `system` coordinates only change over time when `system.ismoving[]`.
-    # `neighbor` coordinates only change over time when `neighbor.ismoving[]`.
+                     u_system, u_neighbor, semi;
+                     # `system` coordinates only change over time when `system.ismoving[]`.
+                     # `neighbor` coordinates only change over time when `neighbor.ismoving[]`.
+                     points_moving=(system.ismoving[], neighbor.ismoving[]))
     update!(neighborhood_search,
             current_coordinates(u_system, system),
             current_coordinates(u_neighbor, neighbor),
-            semi, points_moving=(system.ismoving[], neighbor.ismoving[]))
+            semi; points_moving)
 end
 
 function update_nhs!(neighborhood_search,
                      system::DEMSystem, neighbor::DEMSystem,
-                     u_system, u_neighbor, semi)
-    # Both coordinates change over time
+                     u_system, u_neighbor, semi;
+                     # Both coordinates change over time
+                     points_moving=(true, true))
     update!(neighborhood_search,
             current_coordinates(u_system, system),
             current_coordinates(u_neighbor, neighbor),
-            semi, points_moving=(true, true))
+            semi; points_moving)
 end
 
 function update_nhs!(neighborhood_search,
                      system::DEMSystem, neighbor::BoundaryDEMSystem,
-                     u_system, u_neighbor, semi)
-    # DEM coordinates change over time, the boundary coordinates don't
+                     u_system, u_neighbor, semi;
+                     # DEM coordinates change over time, the boundary coordinates don't
+                     points_moving=(true, false))
     update!(neighborhood_search,
             current_coordinates(u_system, system),
             current_coordinates(u_neighbor, neighbor),
-            semi, points_moving=(true, false))
+            semi; points_moving)
 end
 
 function update_nhs!(neighborhood_search,
                      system::BoundarySPHSystem,
                      neighbor::FluidSystem,
-                     u_system, u_neighbor, semi)
+                     u_system, u_neighbor, semi;
+                     points_moving=(false, false))
     # Don't update. This NHS is never used.
     return neighborhood_search
 end
@@ -814,7 +870,8 @@ end
 function update_nhs!(neighborhood_search,
                      system::BoundaryDEMSystem,
                      neighbor::Union{DEMSystem, BoundaryDEMSystem},
-                     u_system, u_neighbor, semi)
+                     u_system, u_neighbor, semi;
+                     points_moving=(false, false))
     # Don't update. This NHS is never used.
     return neighborhood_search
 end
@@ -822,13 +879,14 @@ end
 function update_nhs!(neighborhood_search,
                      system::Union{BoundarySPHSystem, OpenBoundarySPHSystem},
                      neighbor::Union{BoundarySPHSystem, OpenBoundarySPHSystem},
-                     u_system, u_neighbor, semi)
+                     u_system, u_neighbor, semi;
+                     points_moving=(false, false))
     # Don't update. This NHS is never used.
     return neighborhood_search
 end
 
 # Forward to PointNeighbors.jl
-function update!(neighborhood_search, x, y, semi; points_moving=(true, false))
+function update!(neighborhood_search, x, y, semi; points_moving=(true, true))
     PointNeighbors.update!(neighborhood_search, x, y; points_moving,
                            parallelization_backend=semi.parallelization_backend)
 end
diff --git a/src/schemes/solid/total_lagrangian_sph/rhs.jl b/src/schemes/solid/total_lagrangian_sph/rhs.jl
index 571517799f..9f651aadd6 100644
--- a/src/schemes/solid/total_lagrangian_sph/rhs.jl
+++ b/src/schemes/solid/total_lagrangian_sph/rhs.jl
@@ -2,7 +2,11 @@
 function interact!(dv, v_particle_system, u_particle_system,
                    v_neighbor_system, u_neighbor_system,
                    particle_system::TotalLagrangianSPHSystem,
-                   neighbor_system::TotalLagrangianSPHSystem, semi)
+                   neighbor_system::TotalLagrangianSPHSystem, semi;
+                   integrate_tlsph=semi.integrate_tlsph[])
+    # Skip interaction if TLSPH systems are integrated separately
+    integrate_tlsph || return dv
+
     interact_solid_solid!(dv, particle_system, neighbor_system, semi)
 end
 
@@ -62,7 +66,11 @@ end
 function interact!(dv, v_particle_system, u_particle_system,
                    v_neighbor_system, u_neighbor_system,
                    particle_system::TotalLagrangianSPHSystem,
-                   neighbor_system::FluidSystem, semi)
+                   neighbor_system::FluidSystem, semi;
+                   integrate_tlsph=semi.integrate_tlsph[])
+    # Skip interaction if TLSPH systems are integrated separately
+    integrate_tlsph || return dv
+
     sound_speed = system_sound_speed(neighbor_system)
 
     system_coords = current_coordinates(u_particle_system, particle_system)
@@ -164,7 +172,8 @@ end
 function interact!(dv, v_particle_system, u_particle_system,
                    v_neighbor_system, u_neighbor_system,
                    particle_system::TotalLagrangianSPHSystem,
-                   neighbor_system::Union{BoundarySPHSystem, OpenBoundarySPHSystem}, semi)
+                   neighbor_system::Union{BoundarySPHSystem, OpenBoundarySPHSystem}, semi;
+                   integrate_tlsph=semi.integrate_tlsph[])
     # TODO continuity equation?
     return dv
 end

From 89337b7f1809f8c1e0abcf335f19bddcb0f83ee4 Mon Sep 17 00:00:00 2001
From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com>
Date: Mon, 21 Apr 2025 23:11:01 +0200
Subject: [PATCH 2/9] Update timers

---
 src/callbacks/split_integration.jl | 115 +++++++++++++++--------------
 1 file changed, 60 insertions(+), 55 deletions(-)

diff --git a/src/callbacks/split_integration.jl b/src/callbacks/split_integration.jl
index adee792758..369394a7fb 100644
--- a/src/callbacks/split_integration.jl
+++ b/src/callbacks/split_integration.jl
@@ -65,10 +65,15 @@ function initialize_split_integration!(cb, u, t, integrator)
     ode_split = DynamicalODEProblem(kick_split!, drift_split!, v0_ode_split, u0_ode_split,
                                     tspan, p)
 
-    # Create the split integrator
-    split_integrator = SciMLBase.init(ode_split, alg; save_everystep=false, kwargs...)
+    # Create the split integrator.
+    # We need the timer here to keep the output clean because this will call `kick!` once.
+    @trixi_timeit timer() "split integration" begin
+        split_integrator = SciMLBase.init(ode_split, alg; save_everystep=false, kwargs...)
+    end
+
     # A zero `tspan` sets `tdir` to zero, which breaks adding tstops
     SciMLBase.pop_tstop!(split_integrator)
+
     split_integration_callback.integrator = split_integrator
 
     return cb
@@ -95,7 +100,9 @@ function (split_integration_callback::SplitIntegrationCallback)(integrator)
     @trixi_timeit timer() "split integration" 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.
-        update_systems_and_nhs(v_ode, u_ode, semi, new_t; update_from_callback=true)
+        @trixi_timeit timer() "update systems and nhs" update_systems_and_nhs(v_ode, u_ode,
+                                                                              semi, new_t;
+                                                                              update_from_callback=true)
 
         # Integrate the split integrator to the new time
         add_tstop!(split_integrator, new_t)
@@ -131,68 +138,66 @@ end
 function kick_split!(dv_ode_split, v_ode_split, u_ode_split, p, t)
     (; v_ode, u_ode, semi, semi_split) = p
 
-    @trixi_timeit timer() "split integration kick!" begin
-        @trixi_timeit timer() "reset ∂v/∂t" set_zero!(dv_ode_split)
+    @trixi_timeit timer() "reset ∂v/∂t" set_zero!(dv_ode_split)
 
-        # Only loop over TLSPH systems
-        @trixi_timeit timer() "copy to large v,u" foreach_system(semi_split) do system
-            v = wrap_v(v_ode, system, semi)
-            u = wrap_u(u_ode, system, semi)
-            v_split = wrap_v(v_ode_split, system, semi_split)
-            u_split = wrap_u(u_ode_split, system, semi_split)
+    # Only loop over TLSPH systems
+    @trixi_timeit timer() "copy to large v,u" foreach_system(semi_split) do system
+        v = wrap_v(v_ode, system, semi)
+        u = wrap_u(u_ode, system, semi)
+        v_split = wrap_v(v_ode_split, system, semi_split)
+        u_split = wrap_u(u_ode_split, system, semi_split)
 
-            @threaded semi for particle in each_moving_particle(system)
-                for i in axes(v, 1)
-                    v[i, particle] = v_split[i, particle]
-                end
+        @threaded semi for particle in each_moving_particle(system)
+            for i in axes(v, 1)
+                v[i, particle] = v_split[i, particle]
+            end
 
-                for i in axes(u, 1)
-                    u[i, particle] = u_split[i, particle]
-                end
+            for i in axes(u, 1)
+                u[i, particle] = u_split[i, particle]
             end
         end
+    end
 
-        update_nhs_fun = (semi, u_ode) -> update_nhs_split!(semi, u_ode, u_ode_split,
-                                                            semi_split)
-        @trixi_timeit timer() "update systems and nhs" update_systems_and_nhs(v_ode, u_ode,
-                                                                              semi, t;
-                                                                              systems=semi_split.systems,
-                                                                              update_nhs_fun)
-
-        @trixi_timeit timer() "system interaction" begin
-            # Only loop over systems in the split integrator
-            foreach_system(semi_split) do system
-                # Loop over neighbors in the big integrator
-                foreach_system(semi) do neighbor
-                    # Construct string for the interactions timer.
-                    # Avoid allocations from string construction when no timers are used.
-                    if timeit_debug_enabled()
-                        system_index = system_indices(system, semi)
-                        neighbor_index = system_indices(neighbor, semi)
-                        timer_str = "$(timer_name(system))$system_index-$(timer_name(neighbor))$neighbor_index"
-                    else
-                        timer_str = ""
-                    end
-
-                    dv = wrap_v(dv_ode_split, system, semi_split)
-                    v_system = wrap_v(v_ode_split, system, semi_split)
-                    u_system = wrap_u(u_ode_split, system, semi_split)
-
-                    v_neighbor = wrap_v(v_ode, neighbor, semi)
-                    u_neighbor = wrap_u(u_ode, neighbor, semi)
-
-                    @trixi_timeit timer() timer_str begin
-                        interact!(dv, v_system, u_system, v_neighbor, u_neighbor,
-                                  system, neighbor, semi; integrate_tlsph=true)
-                    end
+    update_nhs_fun = (semi, u_ode) -> update_nhs_split!(semi, u_ode, u_ode_split,
+                                                        semi_split)
+    @trixi_timeit timer() "update systems and nhs" update_systems_and_nhs(v_ode, u_ode,
+                                                                          semi, t;
+                                                                          systems=semi_split.systems,
+                                                                          update_nhs_fun)
+
+    @trixi_timeit timer() "system interaction" begin
+        # Only loop over systems in the split integrator
+        foreach_system(semi_split) do system
+            # Loop over neighbors in the big integrator
+            foreach_system(semi) do neighbor
+                # Construct string for the interactions timer.
+                # Avoid allocations from string construction when no timers are used.
+                if timeit_debug_enabled()
+                    system_index = system_indices(system, semi)
+                    neighbor_index = system_indices(neighbor, semi)
+                    timer_str = "$(timer_name(system))$system_index-$(timer_name(neighbor))$neighbor_index"
+                else
+                    timer_str = ""
+                end
+
+                dv = wrap_v(dv_ode_split, system, semi_split)
+                v_system = wrap_v(v_ode_split, system, semi_split)
+                u_system = wrap_u(u_ode_split, system, semi_split)
+
+                v_neighbor = wrap_v(v_ode, neighbor, semi)
+                u_neighbor = wrap_u(u_ode, neighbor, semi)
+
+                @trixi_timeit timer() timer_str begin
+                    interact!(dv, v_system, u_system, v_neighbor, u_neighbor,
+                              system, neighbor, semi; integrate_tlsph=true)
                 end
             end
         end
-
-        @trixi_timeit timer() "source terms" add_source_terms!(dv_ode_split, v_ode_split,
-                                                               u_ode_split, semi, t;
-                                                               semi_wrap=semi_split)
     end
+
+    @trixi_timeit timer() "source terms" add_source_terms!(dv_ode_split, v_ode_split,
+                                                           u_ode_split, semi, t;
+                                                           semi_wrap=semi_split)
 end
 
 function drift_split!(du_ode, v_ode, u_ode, p, t)

From 524014a66e576a494cb588049594ea7f25e3210e Mon Sep 17 00:00:00 2001
From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com>
Date: Wed, 28 May 2025 11:01:19 +0200
Subject: [PATCH 3/9] Improve split integration callback

---
 src/callbacks/split_integration.jl | 108 +++++++++++++++--------------
 src/general/semidiscretization.jl  |   2 +-
 2 files changed, 57 insertions(+), 53 deletions(-)

diff --git a/src/callbacks/split_integration.jl b/src/callbacks/split_integration.jl
index 369394a7fb..e9f9aac8ac 100644
--- a/src/callbacks/split_integration.jl
+++ b/src/callbacks/split_integration.jl
@@ -42,24 +42,9 @@ function initialize_split_integration!(cb, u, t, integrator)
     u0_ode_split = similar(u_ode, sum(sizes_u))
 
     # Copy the initial conditions to the split integrator
-    foreach_system(semi_split) do system
-        v = wrap_v(v_ode, system, semi)
-        u = wrap_u(u_ode, system, semi)
-        v_split = wrap_v(v0_ode_split, system, semi_split)
-        u_split = wrap_u(u0_ode_split, system, semi_split)
-
-        @threaded semi for particle in each_moving_particle(system)
-            for i in axes(v, 1)
-                v_split[i, particle] = v[i, particle]
-            end
+    copy_to_split!(v_ode, u_ode, v0_ode_split, u0_ode_split, semi, semi_split)
 
-            for i in axes(u, 1)
-                u_split[i, particle] = u[i, particle]
-            end
-        end
-    end
-
-    # No tstops by default
+    # A zero `tspan` sets `tdir` to zero, which breaks adding tstops
     tspan = (integrator.t, integrator.t + 1)
     p = (; v_ode, u_ode, semi, semi_split)
     ode_split = DynamicalODEProblem(kick_split!, drift_split!, v0_ode_split, u0_ode_split,
@@ -71,9 +56,10 @@ function initialize_split_integration!(cb, u, t, integrator)
         split_integrator = SciMLBase.init(ode_split, alg; save_everystep=false, kwargs...)
     end
 
-    # A zero `tspan` sets `tdir` to zero, which breaks adding tstops
+    # Remove the `tstop` (equivalent to zero `tspan`)
     SciMLBase.pop_tstop!(split_integrator)
 
+    # Store the integrator in the callback
     split_integration_callback.integrator = split_integrator
 
     return cb
@@ -81,12 +67,17 @@ end
 
 # `condition`
 function (split_integration_callback::SplitIntegrationCallback)(u, t, integrator)
+    # Integrate TLSPH after every time step
     return true
 end
 
 # `affect!`
 function (split_integration_callback::SplitIntegrationCallback)(integrator)
-    split_integrator = split_integration_callback.integrator
+    # Function barrier for type stability
+    affect_inner!(integrator, split_integration_callback.integrator)
+end
+
+function affect_inner!(integrator, split_integrator)
     semi_split = split_integrator.p.semi_split
 
     semi = integrator.p
@@ -111,22 +102,7 @@ function (split_integration_callback::SplitIntegrationCallback)(integrator)
         v_ode_split, u_ode_split = split_integrator.u.x
 
         # Copy the solutions back to the original integrator
-        foreach_system(semi_split) do system
-            v = wrap_v(v_ode, system, semi)
-            u = wrap_u(u_ode, system, semi)
-            v_split = wrap_v(v_ode_split, system, semi_split)
-            u_split = wrap_u(u_ode_split, system, semi_split)
-
-            @threaded semi for particle in each_moving_particle(system)
-                for i in axes(v, 1)
-                    v[i, particle] = v_split[i, particle]
-                end
-
-                for i in axes(u, 1)
-                    u[i, particle] = u_split[i, particle]
-                end
-            end
-        end
+        copy_from_split!(v_ode, u_ode, v_ode_split, u_ode_split, semi, semi_split)
     end
 
     # Tell OrdinaryDiffEq that `u` has been modified
@@ -141,23 +117,11 @@ function kick_split!(dv_ode_split, v_ode_split, u_ode_split, p, t)
     @trixi_timeit timer() "reset ∂v/∂t" set_zero!(dv_ode_split)
 
     # Only loop over TLSPH systems
-    @trixi_timeit timer() "copy to large v,u" foreach_system(semi_split) do system
-        v = wrap_v(v_ode, system, semi)
-        u = wrap_u(u_ode, system, semi)
-        v_split = wrap_v(v_ode_split, system, semi_split)
-        u_split = wrap_u(u_ode_split, system, semi_split)
-
-        @threaded semi for particle in each_moving_particle(system)
-            for i in axes(v, 1)
-                v[i, particle] = v_split[i, particle]
-            end
-
-            for i in axes(u, 1)
-                u[i, particle] = u_split[i, particle]
-            end
-        end
-    end
+    @trixi_timeit timer() "copy to large v,u" copy_from_split!(v_ode, u_ode,
+                                                               v_ode_split, u_ode_split,
+                                                               semi, semi_split)
 
+    # Update the TLSPH systems with the other systems as neighbors
     update_nhs_fun = (semi, u_ode) -> update_nhs_split!(semi, u_ode, u_ode_split,
                                                         semi_split)
     @trixi_timeit timer() "update systems and nhs" update_systems_and_nhs(v_ode, u_ode,
@@ -168,7 +132,7 @@ function kick_split!(dv_ode_split, v_ode_split, u_ode_split, p, t)
     @trixi_timeit timer() "system interaction" begin
         # Only loop over systems in the split integrator
         foreach_system(semi_split) do system
-            # Loop over neighbors in the big integrator
+            # Loop over all neighbors in the big integrator
             foreach_system(semi) do neighbor
                 # Construct string for the interactions timer.
                 # Avoid allocations from string construction when no timers are used.
@@ -226,6 +190,46 @@ function update_nhs_split!(semi, u_ode, u_ode_split, semi_split)
     end
 end
 
+# Copy the solution from the large integrator to the split integrator
+@inline function copy_to_split!(v_ode, u_ode, v_ode_split, u_ode_split, semi, semi_split)
+    foreach_system(semi_split) do system
+        v = wrap_v(v_ode, system, semi)
+        u = wrap_u(u_ode, system, semi)
+        v_split = wrap_v(v_ode_split, system, semi_split)
+        u_split = wrap_u(u_ode_split, system, semi_split)
+
+        @threaded semi for particle in each_moving_particle(system)
+            for i in axes(v, 1)
+                v_split[i, particle] = v[i, particle]
+            end
+
+            for i in axes(u, 1)
+                u_split[i, particle] = u[i, particle]
+            end
+        end
+    end
+end
+
+# Copy the solution from the split integrator to the large integrator
+@inline function copy_from_split!(v_ode, u_ode, v_ode_split, u_ode_split, semi, semi_split)
+    foreach_system(semi_split) do system
+        v = wrap_v(v_ode, system, semi)
+        u = wrap_u(u_ode, system, semi)
+        v_split = wrap_v(v_ode_split, system, semi_split)
+        u_split = wrap_u(u_ode_split, system, semi_split)
+
+        @threaded semi for particle in each_moving_particle(system)
+            for i in axes(v, 1)
+                v[i, particle] = v_split[i, particle]
+            end
+
+            for i in axes(u, 1)
+                u[i, particle] = u_split[i, particle]
+            end
+        end
+    end
+end
+
 function Base.show(io::IO, cb::DiscreteCallback{<:Any, <:SplitIntegrationCallback})
     @nospecialize cb # reduce precompilation time
     print(io, "SplitIntegrationCallback(integrator=", cb.affect!.integrator, ")")
diff --git a/src/general/semidiscretization.jl b/src/general/semidiscretization.jl
index c587ebbd1d..df835aceeb 100644
--- a/src/general/semidiscretization.jl
+++ b/src/general/semidiscretization.jl
@@ -53,7 +53,7 @@ struct Semidiscretization{BACKEND, S, RU, RV, NS, IT}
     ranges_v                :: RV
     neighborhood_searches   :: NS
     parallelization_backend :: BACKEND
-    integrate_tlsph         :: IT
+    integrate_tlsph         :: IT # `false` if TLSPH integration is decoupled
 
     # Dispatch at `systems` to distinguish this constructor from the one below when
     # 4 systems are passed.

From cd46028599248b9e39effae09e1f212d4f64b14b Mon Sep 17 00:00:00 2001
From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com>
Date: Wed, 28 May 2025 11:05:04 +0200
Subject: [PATCH 4/9] Reformat

---
 src/callbacks/split_integration.jl | 7 ++++---
 src/general/semidiscretization.jl  | 2 +-
 2 files changed, 5 insertions(+), 4 deletions(-)

diff --git a/src/callbacks/split_integration.jl b/src/callbacks/split_integration.jl
index e9f9aac8ac..92fb60c8f8 100644
--- a/src/callbacks/split_integration.jl
+++ b/src/callbacks/split_integration.jl
@@ -14,7 +14,7 @@ function SplitIntegrationCallback(alg; kwargs...)
 
     # The first one is the `condition`, the second the `affect!`
     return DiscreteCallback(split_integration_callback, split_integration_callback,
-                            initialize=initialize_split_integration!,
+                            initialize=(initialize_split_integration!),
                             save_positions=(false, false))
 end
 
@@ -122,8 +122,9 @@ function kick_split!(dv_ode_split, v_ode_split, u_ode_split, p, t)
                                                                semi, semi_split)
 
     # Update the TLSPH systems with the other systems as neighbors
-    update_nhs_fun = (semi, u_ode) -> update_nhs_split!(semi, u_ode, u_ode_split,
-                                                        semi_split)
+    update_nhs_fun = (semi,
+                      u_ode) -> update_nhs_split!(semi, u_ode, u_ode_split,
+                                                  semi_split)
     @trixi_timeit timer() "update systems and nhs" update_systems_and_nhs(v_ode, u_ode,
                                                                           semi, t;
                                                                           systems=semi_split.systems,
diff --git a/src/general/semidiscretization.jl b/src/general/semidiscretization.jl
index df835aceeb..c01bd0559d 100644
--- a/src/general/semidiscretization.jl
+++ b/src/general/semidiscretization.jl
@@ -498,7 +498,7 @@ end
 # before calling `interact!` to compute forces.
 # `semi` and `update_nhs_fun` are overwritten by the `SplitIntegrationCallback`.
 function update_systems_and_nhs(v_ode, u_ode, semi, t;
-                                systems=semi, update_nhs_fun=update_nhs!,
+                                systems=semi, update_nhs_fun=(update_nhs!),
                                 update_from_callback=false)
     # First update step before updating the NHS
     # (for example for writing the current coordinates in the solid system)

From eb897306c6dfb25f7f6ed7a668122580c3f7d214 Mon Sep 17 00:00:00 2001
From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com>
Date: Wed, 28 May 2025 11:35:38 +0200
Subject: [PATCH 5/9] Fix tests

---
 test/test_util.jl | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/test/test_util.jl b/test/test_util.jl
index a487a2d49e..dab317bad8 100644
--- a/test/test_util.jl
+++ b/test/test_util.jl
@@ -45,10 +45,11 @@ macro trixi_testset(name, expr)
 end
 
 struct DummySemidiscretization
-    parallelization_backend::Any
+    parallelization_backend :: Any
+    integrate_tlsph         :: Any
 
     function DummySemidiscretization(; parallelization_backend=SerialBackend())
-        new(parallelization_backend)
+        new(parallelization_backend, Ref(true))
     end
 end
 

From b0e65cbc4a78a65b34dda40f5b6f2e2fdc5951ed Mon Sep 17 00:00:00 2001
From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com>
Date: Wed, 28 May 2025 13:33:03 +0200
Subject: [PATCH 6/9] Fix tests

---
 test/count_allocations.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/test/count_allocations.jl b/test/count_allocations.jl
index eebcaa24ba..0c3bd29533 100644
--- a/test/count_allocations.jl
+++ b/test/count_allocations.jl
@@ -15,7 +15,7 @@ function copy_semi_with_no_update_nhs(semi)
                                   for searches in semi.neighborhood_searches)
 
     return Semidiscretization(semi.systems, semi.ranges_u, semi.ranges_v,
-                              neighborhood_searches, SerialBackend())
+                              neighborhood_searches, SerialBackend(), Ref(true))
 end
 
 # Forward `foreach_neighbor` to wrapped neighborhood search

From f1052026588cb87e8659e757b2c58b46cec693e4 Mon Sep 17 00:00:00 2001
From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com>
Date: Wed, 28 May 2025 14:05:18 +0200
Subject: [PATCH 7/9] Remove a bunch of warnings in the tests

---
 test/schemes/fluid/rhs.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/test/schemes/fluid/rhs.jl b/test/schemes/fluid/rhs.jl
index 21305b73c3..d013d7da33 100644
--- a/test/schemes/fluid/rhs.jl
+++ b/test/schemes/fluid/rhs.jl
@@ -2,6 +2,7 @@
     @testset verbose=true "`pressure_acceleration`" begin
         # Use `@trixi_testset` to isolate the mock functions in a separate namespace
         @trixi_testset "Symmetry" begin
+            TrixiParticles.ndims(::Val{:smoothing_kernel}) = 2
             masses = [[0.01, 0.01], [0.73, 0.31]]
             densities = [
                 [1000.0, 1000.0],
@@ -47,7 +48,6 @@
                         density = ones(3)
                         state_equation = Val(:state_equation)
                         smoothing_kernel = Val(:smoothing_kernel)
-                        TrixiParticles.ndims(::Val{:smoothing_kernel}) = 2
                         smoothing_length = -1.0
 
                         fluid = InitialCondition(; coordinates, velocity, mass, density)

From 794cbfb13bfd64b7de1825f609de6b98e9d7c026 Mon Sep 17 00:00:00 2001
From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com>
Date: Wed, 28 May 2025 18:08:48 +0200
Subject: [PATCH 8/9] Improve timers

---
 src/TrixiParticles.jl              |  2 +-
 src/callbacks/split_integration.jl | 12 +++++++-----
 2 files changed, 8 insertions(+), 6 deletions(-)

diff --git a/src/TrixiParticles.jl b/src/TrixiParticles.jl
index e90574aeb7..eeec9a02b4 100644
--- a/src/TrixiParticles.jl
+++ b/src/TrixiParticles.jl
@@ -28,7 +28,7 @@ using SciMLBase: SciMLBase, CallbackSet, DiscreteCallback, DynamicalODEProblem,
 @reexport using StaticArrays: SVector
 using StaticArrays: @SMatrix, SMatrix, setindex
 using StrideArrays: PtrArray, StaticInt
-using TimerOutputs: TimerOutput, TimerOutputs, print_timer, reset_timer!
+using TimerOutputs: TimerOutput, TimerOutputs, print_timer, reset_timer!, @notimeit
 using TrixiBase: @trixi_timeit, timer, timeit_debug_enabled,
                  disable_debug_timings, enable_debug_timings, TrixiBase
 @reexport using TrixiBase: trixi_include, trixi_include_changeprecision
diff --git a/src/callbacks/split_integration.jl b/src/callbacks/split_integration.jl
index 92fb60c8f8..441fcf0290 100644
--- a/src/callbacks/split_integration.jl
+++ b/src/callbacks/split_integration.jl
@@ -53,7 +53,9 @@ function initialize_split_integration!(cb, u, t, integrator)
     # Create the split integrator.
     # We need the timer here to keep the output clean because this will call `kick!` once.
     @trixi_timeit timer() "split integration" begin
-        split_integrator = SciMLBase.init(ode_split, alg; save_everystep=false, kwargs...)
+        @trixi_timeit timer() "init" begin
+            TimerOutputs.@notimeit timer() split_integrator = SciMLBase.init(ode_split, alg; save_everystep=false, kwargs...)
+        end
     end
 
     # Remove the `tstop` (equivalent to zero `tspan`)
@@ -97,12 +99,12 @@ function affect_inner!(integrator, split_integrator)
 
         # Integrate the split integrator to the new time
         add_tstop!(split_integrator, new_t)
-        SciMLBase.solve!(split_integrator)
+        @trixi_timeit timer() "solve" SciMLBase.solve!(split_integrator)
 
         v_ode_split, u_ode_split = split_integrator.u.x
 
         # Copy the solutions back to the original integrator
-        copy_from_split!(v_ode, u_ode, v_ode_split, u_ode_split, semi, semi_split)
+        @trixi_timeit timer() "copy back" copy_from_split!(v_ode, u_ode, v_ode_split, u_ode_split, semi, semi_split)
     end
 
     # Tell OrdinaryDiffEq that `u` has been modified
@@ -233,7 +235,7 @@ end
 
 function Base.show(io::IO, cb::DiscreteCallback{<:Any, <:SplitIntegrationCallback})
     @nospecialize cb # reduce precompilation time
-    print(io, "SplitIntegrationCallback(integrator=", cb.affect!.integrator, ")")
+    print(io, "SplitIntegrationCallback(alg=", cb.affect!.alg, ")")
 end
 
 function Base.show(io::IO, ::MIME"text/plain",
@@ -245,7 +247,7 @@ function Base.show(io::IO, ::MIME"text/plain",
     else
         update_cb = cb.affect!
         setup = [
-            "integrator" => update_cb.integrator
+            "alg" => update_cb.alg
         ]
         summary_box(io, "SplitIntegrationCallback", setup)
     end

From da28e75ca29e9fbd8d370d26fbd3e51fe617f244 Mon Sep 17 00:00:00 2001
From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com>
Date: Mon, 23 Jun 2025 15:21:40 +0200
Subject: [PATCH 9/9] Fix split integration on GPUs

---
 src/general/semidiscretization.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/general/semidiscretization.jl b/src/general/semidiscretization.jl
index 57100a638f..b5ffe0c632 100644
--- a/src/general/semidiscretization.jl
+++ b/src/general/semidiscretization.jl
@@ -330,7 +330,7 @@ function semidiscretize(semi, tspan; reset_threads=true)
         semi_new = Semidiscretization(set_system_links.(semi_.systems, Ref(semi_)),
                                       semi_.ranges_u, semi_.ranges_v,
                                       semi_.neighborhood_searches,
-                                      semi_.parallelization_backend)
+                                      semi_.parallelization_backend, semi_.integrate_tlsph)
     else
         semi_new = semi
     end