From 9b01feef60fd7ecb1843cff5584c3ecb172074bb Mon Sep 17 00:00:00 2001
From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com>
Date: Wed, 24 Sep 2025 13:20:43 +0200
Subject: [PATCH 01/16] Add `EnergyCalculatorCallback`

---
 src/callbacks/callbacks.jl                    |   1 +
 src/callbacks/energy_calculator.jl            | 158 ++++++++++++++++++
 src/callbacks/update.jl                       |   2 +-
 .../structure/total_lagrangian_sph/rhs.jl     |  30 ++--
 .../structure/total_lagrangian_sph/system.jl  |  14 +-
 5 files changed, 189 insertions(+), 16 deletions(-)
 create mode 100644 src/callbacks/energy_calculator.jl

diff --git a/src/callbacks/callbacks.jl b/src/callbacks/callbacks.jl
index b1bd602269..1be909f86a 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("split_integration.jl")
+include("energy_calculator.jl")
diff --git a/src/callbacks/energy_calculator.jl b/src/callbacks/energy_calculator.jl
new file mode 100644
index 0000000000..be1ebdef58
--- /dev/null
+++ b/src/callbacks/energy_calculator.jl
@@ -0,0 +1,158 @@
+struct EnergyCalculatorCallback{T}
+    interval :: Int
+    t        :: T
+    energy   :: T
+end
+
+function EnergyCalculatorCallback{ELTYPE}(; interval=1) where {ELTYPE <: Real}
+    cb = EnergyCalculatorCallback(interval, Ref(zero(ELTYPE)), Ref(zero(ELTYPE)))
+
+    # The first one is the `condition`, the second the `affect!`
+    return DiscreteCallback(cb, cb, save_positions=(false, false))
+end
+
+# `condition`
+function (callback::EnergyCalculatorCallback)(u, t, integrator)
+    (; interval) = callback
+
+    return condition_integrator_interval(integrator, interval)
+end
+
+# `affect!`
+function (callback::EnergyCalculatorCallback)(integrator)
+    # Determine time step size as difference to last time this callback was called
+    t = integrator.t
+    dt = t - callback.t[]
+
+    # Update time of last call
+    callback.t[] = t
+
+    semi = integrator.p
+    v_ode, u_ode = integrator.u.x
+    energy = callback.energy
+
+    update_energy_calculator!(energy, v_ode, u_ode, semi.systems[1], semi, t, dt)
+
+    # Tell OrdinaryDiffEq that `u` has not been modified
+    u_modified!(integrator, false)
+
+    return integrator
+end
+
+function update_energy_calculator!(energy, v_ode, u_ode,
+                                   system::AbstractStructureSystem, semi, t, dt)
+    @trixi_timeit timer() "calculate energy" begin
+        # Update quantities that are stored in the systems. These quantities (e.g. pressure)
+        # still have the values from the last stage of the previous step if not updated here.
+        @trixi_timeit timer() "update systems and nhs" begin
+            # Don't create sub-timers here to avoid cluttering the timer output
+            @notimeit timer() update_systems_and_nhs(v_ode, u_ode, semi, t)
+        end
+
+        dv_clamped = system.cache.dv_clamped
+        set_zero!(dv_clamped)
+
+        dv = zero(wrap_v(v_ode, system, semi))
+        dv_combined = CombinedMatrix(dv, dv_clamped)
+        v = wrap_v(v_ode, system, semi)
+        u = wrap_u(u_ode, system, semi)
+        eachparticle = (n_integrated_particles(system) + 1):nparticles(system)
+
+        foreach_system(semi) do neighbor_system
+            v_neighbor = wrap_v(v_ode, neighbor_system, semi)
+            u_neighbor = wrap_u(u_ode, neighbor_system, semi)
+
+            interact!(dv_combined, v, u, v_neighbor, u_neighbor,
+                      system, neighbor_system, semi,
+                      eachparticle=eachparticle)
+        end
+
+        @threaded semi for particle in eachparticle
+            # Dispatch by system type to exclude boundary systems
+            add_acceleration!(dv_combined, particle, system)
+            add_source_terms_inner!(dv_combined, v, u, particle, system,
+                                    source_terms(system), t)
+        end
+
+        n_clamped_particles = nparticles(system) - n_integrated_particles(system)
+        for clamped_particle in 1:n_clamped_particles
+            particle = clamped_particle + n_integrated_particles(system)
+            velocity = current_velocity(nothing, system, particle)
+            dv_particle = extract_svector(dv_clamped, system, clamped_particle)
+
+            # The force on the clamped particle is mass times acceleration
+            F_particle = system.mass[particle] * dv_particle
+
+            # To obtain energy, we need to integrate the instantaneous power.
+            # Instantaneous power is force done BY the particle times prescribed velocity.
+            # The work done BY the particle is the negative of the work done ON it.
+            energy[] -= dot(F_particle, velocity) * dt
+        end
+    end
+end
+
+# Data type that combined two matrices and behaves like a single matrix.
+# This is used in `TotalLagrangianSPHSystem` to combine the `dv` for moving
+# particles and the `dv_fixed` for fixed particles into a single matrix that
+# can be passed to `interact!`.
+# This is required when `compute_forces_for_fixed_particles=true` is used.
+struct CombinedMatrix{T, M1, M2} <: AbstractMatrix{T}
+    matrix1::M1
+    matrix2::M2
+
+    function CombinedMatrix(matrix1, matrix2)
+        @assert size(matrix1, 1) == size(matrix2, 1)
+        @assert eltype(matrix1) == eltype(matrix2)
+
+        new{eltype(matrix1), typeof(matrix1), typeof(matrix2)}(matrix1, matrix2)
+    end
+end
+
+@inline function Base.size(cm::CombinedMatrix)
+    return (size(cm.matrix1, 1), size(cm.matrix1, 2) + size(cm.matrix2, 2))
+end
+
+@inline function Base.getindex(cm::CombinedMatrix, i, j)
+    @boundscheck checkbounds(cm, i, j)
+
+    length1 = size(cm.matrix1, 2)
+    if j <= length1
+        return @inbounds cm.matrix1[i, j]
+    else
+        return @inbounds cm.matrix2[i, j - length1]
+    end
+end
+
+@inline function Base.setindex!(cm::CombinedMatrix, value, i, j)
+    @boundscheck checkbounds(cm, i, j)
+
+    length1 = size(cm.matrix1, 2)
+    if j <= length1
+        return @inbounds cm.matrix1[i, j] = value
+    else
+        return @inbounds cm.matrix2[i, j - length1] = value
+    end
+end
+
+function Base.show(io::IO, cb::DiscreteCallback{<:Any, <:EnergyCalculatorCallback})
+    @nospecialize cb # reduce precompilation time
+
+    ELTYPE = eltype(cb.affect!.energy)
+    print(io, "EnergyCalculatorCallback{$ELTYPE}(interval=", cb.affect!.interval, ")")
+end
+
+function Base.show(io::IO, ::MIME"text/plain",
+                   cb::DiscreteCallback{<:Any, EnergyCalculatorCallback{T}}) where {T}
+    @nospecialize cb # reduce precompilation time
+
+    if get(io, :compact, false)
+        show(io, cb)
+    else
+        update_cb = cb.affect!
+        ELTYPE = eltype(update_cb.energy)
+        setup = [
+            "interval" => update_cb.interval
+        ]
+        summary_box(io, "EnergyCalculatorCallback{$ELTYPE}", setup)
+    end
+end
diff --git a/src/callbacks/update.jl b/src/callbacks/update.jl
index 8eecb8380a..056191330e 100644
--- a/src/callbacks/update.jl
+++ b/src/callbacks/update.jl
@@ -31,7 +31,7 @@ function UpdateCallback(; interval::Integer=-1, dt=0.0)
     update_callback! = UpdateCallback(interval)
 
     if dt > 0
-        # Add a `tstop` every `dt`, and save the final solution.
+        # Add a `tstop` every `dt`
         return PeriodicCallback(update_callback!, dt,
                                 initialize=(initial_update!),
                                 save_positions=(false, false))
diff --git a/src/schemes/structure/total_lagrangian_sph/rhs.jl b/src/schemes/structure/total_lagrangian_sph/rhs.jl
index 3edad47a55..683bc0bd80 100644
--- a/src/schemes/structure/total_lagrangian_sph/rhs.jl
+++ b/src/schemes/structure/total_lagrangian_sph/rhs.jl
@@ -3,18 +3,21 @@ function interact!(dv, v_particle_system, u_particle_system,
                    v_neighbor_system, u_neighbor_system,
                    particle_system::TotalLagrangianSPHSystem,
                    neighbor_system::TotalLagrangianSPHSystem, semi;
-                   integrate_tlsph=semi.integrate_tlsph[])
+                   integrate_tlsph=semi.integrate_tlsph[],
+                   eachparticle=each_integrated_particle(particle_system))
     # Different structures do not interact with each other (yet)
     particle_system === neighbor_system || return dv
 
     # Skip interaction if TLSPH systems are integrated separately
     integrate_tlsph || return dv
 
-    interact_structure_structure!(dv, v_particle_system, particle_system, semi)
+    interact_structure_structure!(dv, v_particle_system, particle_system, semi;
+                                  eachparticle)
 end
 
 # Function barrier without dispatch for unit testing
-@inline function interact_structure_structure!(dv, v_system, system, semi)
+@inline function interact_structure_structure!(dv, v_system, system, semi;
+                                               eachparticle=each_integrated_particle(system))
     (; penalty_force) = system
 
     # Everything here is done in the initial coordinates
@@ -23,10 +26,10 @@ end
     # Loop over all pairs of particles and neighbors within the kernel cutoff.
     # For structure-structure interaction, this has to happen in the initial coordinates.
     foreach_point_neighbor(system, system, system_coords, system_coords, semi;
-                           points=each_integrated_particle(system)) do particle, neighbor,
-                                                                       initial_pos_diff,
-                                                                       initial_distance
-        # Only consider particles with a distance > 0. See `src/general/smoothing_kernels.jl` for more details.
+                           points=eachparticle) do particle, neighbor,
+                                                   initial_pos_diff, initial_distance
+        # Only consider particles with a distance > 0.
+        # See `src/general/smoothing_kernels.jl` for more details.
         initial_distance^2 < eps(initial_smoothing_length(system)^2) && return
 
         rho_a = @inbounds system.material_density[particle]
@@ -72,7 +75,8 @@ function interact!(dv, v_particle_system, u_particle_system,
                    v_neighbor_system, u_neighbor_system,
                    particle_system::TotalLagrangianSPHSystem,
                    neighbor_system::AbstractFluidSystem, semi;
-                   integrate_tlsph=semi.integrate_tlsph[])
+                   integrate_tlsph=semi.integrate_tlsph[],
+                   eachparticle=each_integrated_particle(particle_system))
     # Skip interaction if TLSPH systems are integrated separately
     integrate_tlsph || return dv
 
@@ -83,11 +87,8 @@ function interact!(dv, v_particle_system, u_particle_system,
 
     # Loop over all pairs of particles and neighbors within the kernel cutoff
     foreach_point_neighbor(particle_system, neighbor_system, system_coords, neighbor_coords,
-                           semi;
-                           points=each_integrated_particle(particle_system)) do particle,
-                                                                                neighbor,
-                                                                                pos_diff,
-                                                                                distance
+                           semi; points=eachparticle) do particle, neighbor,
+                                                         pos_diff, distance
         # Only consider particles with a distance > 0. See `src/general/smoothing_kernels.jl` for more details.
         distance^2 < eps(initial_smoothing_length(particle_system)^2) && return
 
@@ -178,7 +179,8 @@ function interact!(dv, v_particle_system, u_particle_system,
                    v_neighbor_system, u_neighbor_system,
                    particle_system::TotalLagrangianSPHSystem,
                    neighbor_system::Union{WallBoundarySystem, OpenBoundarySystem}, semi;
-                   integrate_tlsph=semi.integrate_tlsph[])
+                   integrate_tlsph=semi.integrate_tlsph[],
+                   eachparticle=each_integrated_particle(particle_system))
     # TODO continuity equation?
     return dv
 end
diff --git a/src/schemes/structure/total_lagrangian_sph/system.jl b/src/schemes/structure/total_lagrangian_sph/system.jl
index 3979093ea0..a8e9d97323 100644
--- a/src/schemes/structure/total_lagrangian_sph/system.jl
+++ b/src/schemes/structure/total_lagrangian_sph/system.jl
@@ -91,6 +91,7 @@ end
 function TotalLagrangianSPHSystem(initial_condition, smoothing_kernel, smoothing_length,
                                   young_modulus, poisson_ratio;
                                   n_clamped_particles=0, clamped_particles_motion=nothing,
+                                  use_with_energy_calculator_callback=false,
                                   boundary_model=nothing,
                                   acceleration=ntuple(_ -> 0.0,
                                                       ndims(smoothing_kernel)),
@@ -128,7 +129,9 @@ function TotalLagrangianSPHSystem(initial_condition, smoothing_kernel, smoothing
     initialize_prescribed_motion!(clamped_particles_motion, initial_condition,
                                   n_clamped_particles)
 
-    cache = create_cache_tlsph(clamped_particles_motion, initial_condition)
+    cache = (; create_cache_tlsph(clamped_particles_motion, initial_condition)...,
+             create_cache_tlsph(Val(use_with_energy_calculator_callback),
+                                initial_condition, n_clamped_particles)...)
 
     return TotalLagrangianSPHSystem(initial_condition, initial_coordinates,
                                     current_coordinates, mass, correction_matrix,
@@ -149,6 +152,15 @@ function create_cache_tlsph(::PrescribedMotion, initial_condition)
     return (; velocity, acceleration)
 end
 
+create_cache_tlsph(::Val{false}, initial_condition, n_clamped_particles) = (;)
+
+function create_cache_tlsph(::Val{true}, initial_condition, n_clamped_particles)
+    dv_clamped = Array{eltype(initial_condition)}(undef, ndims(initial_condition),
+                                                n_clamped_particles)
+
+    return (; dv_clamped)
+end
+
 @inline function Base.eltype(::TotalLagrangianSPHSystem{<:Any, <:Any, ELTYPE}) where {ELTYPE}
     return ELTYPE
 end

From 79eb1f1b05bd180e192a1c7a0160348dad62941d Mon Sep 17 00:00:00 2001
From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com>
Date: Mon, 6 Oct 2025 10:06:14 +0200
Subject: [PATCH 02/16] Improve comments

---
 src/callbacks/energy_calculator.jl | 54 ++++++++++++++----------------
 src/callbacks/split_integration.jl |  7 ++--
 2 files changed, 30 insertions(+), 31 deletions(-)

diff --git a/src/callbacks/energy_calculator.jl b/src/callbacks/energy_calculator.jl
index be1ebdef58..6c47899620 100644
--- a/src/callbacks/energy_calculator.jl
+++ b/src/callbacks/energy_calculator.jl
@@ -91,11 +91,32 @@ function update_energy_calculator!(energy, v_ode, u_ode,
     end
 end
 
-# Data type that combined two matrices and behaves like a single matrix.
-# This is used in `TotalLagrangianSPHSystem` to combine the `dv` for moving
-# particles and the `dv_fixed` for fixed particles into a single matrix that
-# can be passed to `interact!`.
-# This is required when `compute_forces_for_fixed_particles=true` is used.
+function Base.show(io::IO, cb::DiscreteCallback{<:Any, <:EnergyCalculatorCallback})
+    @nospecialize cb # reduce precompilation time
+
+    ELTYPE = eltype(cb.affect!.energy)
+    print(io, "EnergyCalculatorCallback{$ELTYPE}(interval=", cb.affect!.interval, ")")
+end
+
+function Base.show(io::IO, ::MIME"text/plain",
+                   cb::DiscreteCallback{<:Any, EnergyCalculatorCallback{T}}) where {T}
+    @nospecialize cb # reduce precompilation time
+
+    if get(io, :compact, false)
+        show(io, cb)
+    else
+        update_cb = cb.affect!
+        ELTYPE = eltype(update_cb.energy)
+        setup = [
+            "interval" => update_cb.interval
+        ]
+        summary_box(io, "EnergyCalculatorCallback{$ELTYPE}", setup)
+    end
+end
+
+# Data type that combines two matrices to behave like a single matrix.
+# This is used above to combine the `dv` for integrated particles and the `dv_clamped`
+# for clamped particles into a single matrix that can be passed to `interact!`.
 struct CombinedMatrix{T, M1, M2} <: AbstractMatrix{T}
     matrix1::M1
     matrix2::M2
@@ -133,26 +154,3 @@ end
         return @inbounds cm.matrix2[i, j - length1] = value
     end
 end
-
-function Base.show(io::IO, cb::DiscreteCallback{<:Any, <:EnergyCalculatorCallback})
-    @nospecialize cb # reduce precompilation time
-
-    ELTYPE = eltype(cb.affect!.energy)
-    print(io, "EnergyCalculatorCallback{$ELTYPE}(interval=", cb.affect!.interval, ")")
-end
-
-function Base.show(io::IO, ::MIME"text/plain",
-                   cb::DiscreteCallback{<:Any, EnergyCalculatorCallback{T}}) where {T}
-    @nospecialize cb # reduce precompilation time
-
-    if get(io, :compact, false)
-        show(io, cb)
-    else
-        update_cb = cb.affect!
-        ELTYPE = eltype(update_cb.energy)
-        setup = [
-            "interval" => update_cb.interval
-        ]
-        summary_box(io, "EnergyCalculatorCallback{$ELTYPE}", setup)
-    end
-end
diff --git a/src/callbacks/split_integration.jl b/src/callbacks/split_integration.jl
index 3e7b5fce32..d09b090bd9 100644
--- a/src/callbacks/split_integration.jl
+++ b/src/callbacks/split_integration.jl
@@ -1,6 +1,6 @@
 mutable struct SplitIntegrationCallback
     # Type-stability does not matter here, and it is impossible to implement because
-    # the integration will be created during the initialization.
+    # the integrator will be created during the initialization.
     integrator :: Any
     alg        :: Any
     kwargs     :: Any
@@ -71,9 +71,10 @@ function initialize_split_integration!(cb, u, t, integrator)
     # Create split integrator with TLSPH systems only
     systems = filter(i -> i isa TotalLagrangianSPHSystem, semi.systems)
 
-    # These neighborhood searches are never used
+    # This neighborhood search is never used
+    neighborhood_search = TrivialNeighborhoodSearch{ndims(first(systems))}()
     semi_split = Semidiscretization(systems...,
-                                    neighborhood_search=TrivialNeighborhoodSearch{ndims(first(systems))}(),
+                                    neighborhood_search=neighborhood_search,
                                     parallelization_backend=semi.parallelization_backend)
 
     # Verify that a NHS implementation is used that does not require updates

From d15292775f157bb130d689079aaa0d5c1dcc14ce Mon Sep 17 00:00:00 2001
From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com>
Date: Mon, 6 Oct 2025 11:20:39 +0200
Subject: [PATCH 03/16] Add docstring

---
 src/TrixiParticles.jl              |  2 +-
 src/callbacks/energy_calculator.jl | 53 +++++++++++++++++++++++++++++-
 2 files changed, 53 insertions(+), 2 deletions(-)

diff --git a/src/TrixiParticles.jl b/src/TrixiParticles.jl
index bb9160c848..70b733724a 100644
--- a/src/TrixiParticles.jl
+++ b/src/TrixiParticles.jl
@@ -70,7 +70,7 @@ export WeaklyCompressibleSPHSystem, EntropicallyDampedSPHSystem, TotalLagrangian
 export BoundaryZone, InFlow, OutFlow, BidirectionalFlow
 export InfoCallback, SolutionSavingCallback, DensityReinitializationCallback,
        PostprocessCallback, StepsizeCallback, UpdateCallback, SteadyStateReachedCallback,
-       SplitIntegrationCallback
+       SplitIntegrationCallback, EnergyCalculatorCallback
 export ContinuityDensity, SummationDensity
 export PenaltyForceGanzenmueller, TransportVelocityAdami, ParticleShiftingTechnique,
        ParticleShiftingTechniqueSun2017, ConsistentShiftingSun2019,
diff --git a/src/callbacks/energy_calculator.jl b/src/callbacks/energy_calculator.jl
index 6c47899620..7916b643c7 100644
--- a/src/callbacks/energy_calculator.jl
+++ b/src/callbacks/energy_calculator.jl
@@ -1,3 +1,55 @@
+
+"""
+    EnergyCalculatorCallback{ELTYPE}(; interval=1)
+
+Callback to calculate the energy contribution from clamped particles in
+[`TotalLagrangianSPHSystem`](@ref)s.
+
+This callback computes the work done by clamped particles that are moving according to a
+[`PrescribedMotion`](@ref).
+The energy is calculated by integrating the instantaneous power, which is the dot product
+of the force exerted by the clamped particle and its prescribed velocity.
+The callback uses a simple explicit Euler time integration for the energy calculation.
+
+!!! info "TLSPH System Configuration"
+    The [`TotalLagrangianSPHSystem`](@ref) must be created with
+    `use_with_energy_calculator_callback=true` to prepare the system for energy
+    calculation by allocating the necessary cache for clamped particles.
+
+!!! warning "Experimental implementation"
+    This is an experimental feature and may change in future releases.
+
+# Arguments
+- `ELTYPE`: The floating point type used for time and energy calculations.
+            This should match the floating point type used in the
+            [`TotalLagrangianSPHSystem`](@ref), which can be obtained via `eltype(system)`.
+
+# Keywords
+- `interval=1`: Interval (in number of time steps) at which to compute the energy.
+                It is recommended to set this either to `1` (every time step) or to a small
+                integer (e.g., `2` or `5`) to reduce time integration errors
+                in the energy calculation.
+
+# Examples
+```jldoctest; output = false, setup = :(structure = RectangularShape(0.1, (3, 4), (0.1, 0.0), density=1.0); smoothing_kernel = WendlandC2Kernel{2}(); smoothing_length = 1.0; young_modulus = 1e6; poisson_ratio = 0.3; n_clamped_particles = 0; clamped_particles_motion = nothing)
+# Create TLSPH system with `use_with_energy_calculator_callback=true`
+system = TotalLagrangianSPHSystem(structure, smoothing_kernel, smoothing_length,
+                                  young_modulus, poisson_ratio,
+                                  n_clamped_particles=n_clamped_particles,
+                                  clamped_particles_motion=clamped_particles_motion,
+                                  use_with_energy_calculator_callback=true) # Important!
+
+# Create an energy calculator that runs every 2 time steps
+energy_cb = EnergyCalculatorCallback{eltype(system)}(; interval=2)
+
+# output
+┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
+│ EnergyCalculatorCallback{Float64}                                                                │
+│ ═════════════════════════════════                                                                │
+│ interval: ……………………………………………………… 2                                                                │
+└──────────────────────────────────────────────────────────────────────────────────────────────────┘
+```
+"""
 struct EnergyCalculatorCallback{T}
     interval :: Int
     t        :: T
@@ -68,7 +120,6 @@ function update_energy_calculator!(energy, v_ode, u_ode,
         end
 
         @threaded semi for particle in eachparticle
-            # Dispatch by system type to exclude boundary systems
             add_acceleration!(dv_combined, particle, system)
             add_source_terms_inner!(dv_combined, v, u, particle, system,
                                     source_terms(system), t)

From be1b585e03be40ad8efced0a472c584ebd3aa4fb Mon Sep 17 00:00:00 2001
From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com>
Date: Mon, 6 Oct 2025 11:22:25 +0200
Subject: [PATCH 04/16] Reformat

---
 src/schemes/structure/total_lagrangian_sph/rhs.jl    | 5 +++--
 src/schemes/structure/total_lagrangian_sph/system.jl | 2 +-
 2 files changed, 4 insertions(+), 3 deletions(-)

diff --git a/src/schemes/structure/total_lagrangian_sph/rhs.jl b/src/schemes/structure/total_lagrangian_sph/rhs.jl
index 683bc0bd80..498cc78e59 100644
--- a/src/schemes/structure/total_lagrangian_sph/rhs.jl
+++ b/src/schemes/structure/total_lagrangian_sph/rhs.jl
@@ -87,8 +87,9 @@ function interact!(dv, v_particle_system, u_particle_system,
 
     # Loop over all pairs of particles and neighbors within the kernel cutoff
     foreach_point_neighbor(particle_system, neighbor_system, system_coords, neighbor_coords,
-                           semi; points=eachparticle) do particle, neighbor,
-                                                         pos_diff, distance
+                           semi;
+                           points=eachparticle) do particle, neighbor,
+                                                   pos_diff, distance
         # Only consider particles with a distance > 0. See `src/general/smoothing_kernels.jl` for more details.
         distance^2 < eps(initial_smoothing_length(particle_system)^2) && return
 
diff --git a/src/schemes/structure/total_lagrangian_sph/system.jl b/src/schemes/structure/total_lagrangian_sph/system.jl
index a8e9d97323..c9982855f8 100644
--- a/src/schemes/structure/total_lagrangian_sph/system.jl
+++ b/src/schemes/structure/total_lagrangian_sph/system.jl
@@ -156,7 +156,7 @@ create_cache_tlsph(::Val{false}, initial_condition, n_clamped_particles) = (;)
 
 function create_cache_tlsph(::Val{true}, initial_condition, n_clamped_particles)
     dv_clamped = Array{eltype(initial_condition)}(undef, ndims(initial_condition),
-                                                n_clamped_particles)
+                                                  n_clamped_particles)
 
     return (; dv_clamped)
 end

From d54ea66afdf2973ddf3dbb92602d47be713b105b Mon Sep 17 00:00:00 2001
From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com>
Date: Mon, 6 Oct 2025 14:44:15 +0200
Subject: [PATCH 05/16] Add tests

---
 .../structure/total_lagrangian_sph/system.jl  |  10 +-
 test/callbacks/callbacks.jl                   |   1 +
 test/callbacks/energy_calculator.jl           | 124 ++++++++++++++++++
 test/examples/examples.jl                     |  36 +++++
 4 files changed, 168 insertions(+), 3 deletions(-)
 create mode 100644 test/callbacks/energy_calculator.jl

diff --git a/src/schemes/structure/total_lagrangian_sph/system.jl b/src/schemes/structure/total_lagrangian_sph/system.jl
index c9982855f8..c88bd42534 100644
--- a/src/schemes/structure/total_lagrangian_sph/system.jl
+++ b/src/schemes/structure/total_lagrangian_sph/system.jl
@@ -348,11 +348,15 @@ end
     # Reset deformation gradient
     set_zero!(deformation_grad)
 
-    # Loop over all pairs of particles and neighbors within the kernel cutoff.
+    # Loop over all pairs of particles and neighbors within the kernel cutoff
     initial_coords = initial_coordinates(system)
     foreach_point_neighbor(system, system, initial_coords, initial_coords,
-                           semi) do particle, neighbor, initial_pos_diff, initial_distance
-        # Only consider particles with a distance > 0. See `src/general/smoothing_kernels.jl` for more details.
+                           semi,
+                           parallelization_backend=SerialBackend()) do particle, neighbor,
+                                                                       initial_pos_diff,
+                                                                       initial_distance
+        # Only consider particles with a distance > 0.
+        # See `src/general/smoothing_kernels.jl` for more details.
         initial_distance^2 < eps(initial_smoothing_length(system)^2) && return
 
         volume = mass[neighbor] / material_density[neighbor]
diff --git a/test/callbacks/callbacks.jl b/test/callbacks/callbacks.jl
index ae6166cde8..9c1ea952e7 100644
--- a/test/callbacks/callbacks.jl
+++ b/test/callbacks/callbacks.jl
@@ -5,4 +5,5 @@
     include("update.jl")
     include("solution_saving.jl")
     include("steady_state_reached.jl")
+    include("energy_calculator.jl")
 end
diff --git a/test/callbacks/energy_calculator.jl b/test/callbacks/energy_calculator.jl
new file mode 100644
index 0000000000..214efd9bbd
--- /dev/null
+++ b/test/callbacks/energy_calculator.jl
@@ -0,0 +1,124 @@
+@testset verbose=true "EnergyCalculatorCallback" begin
+    @testset "Constructor and Basic Properties" begin
+        # Test default constructor
+        callback = EnergyCalculatorCallback{Float64}()
+        @test callback.affect!.interval == 1
+        @test callback.affect!.t[] == 0.0
+        @test callback.affect!.energy[] == 0.0
+
+        # Test constructor with interval
+        callback = EnergyCalculatorCallback{Float64}(; interval=5)
+        @test callback.affect!.interval == 5
+
+        # Test with specific element type
+        callback = EnergyCalculatorCallback{Float32}(; interval=2)
+        @test eltype(callback.affect!.energy) == Float32
+        @test eltype(callback.affect!.t) == Float32
+    end
+
+    @testset "show" begin
+        callback = EnergyCalculatorCallback{Float64}(; interval=10)
+
+        # Test compact representation
+        show_compact = "EnergyCalculatorCallback{Float64}(interval=10)"
+        @test repr(callback) == show_compact
+
+        # Test detailed representation - check against expected box format
+        show_box = """
+        ┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
+        │ EnergyCalculatorCallback{Float64}                                                                │
+        │ ═════════════════════════════════                                                                │
+        │ interval: ……………………………………………………… 10                                                               │
+        └──────────────────────────────────────────────────────────────────────────────────────────────────┘"""
+        @test repr("text/plain", callback) == show_box
+    end
+
+    @testset "update_energy_calculator!" begin
+        # In the first test, we just move the 2x2 grid of particles up against gravity
+        # and test that the energy calculated is just the potential energy difference.
+        # In the other tests, we clamp the top row of particles and offset them to create
+        # stress. We can then test how much energy is required to pull the particles further
+        # apart against the elastic forces.
+        n_clamped_particles = [4, 2, 2, 2]
+        E = [1e6, 1e1, 1e6, 1e8]
+        @testset "n_clamped_particles = $(n_clamped_particles[i]), E = $(E[i])" for i in
+                                                                                    eachindex(E)
+            # Create a simple 2D system with 2x2 particles
+            coordinates = [0.0 1.0 0.0 1.0; 0.0 0.0 1.0 1.0]
+            mass = [1.0, 1.0, 1.0, 1.0]
+            density = [1000.0, 1000.0, 1000.0, 1000.0]
+
+            smoothing_kernel = SchoenbergCubicSplineKernel{2}()
+            smoothing_length = sqrt(2)
+            young_modulus = E[i]
+            poisson_ratio = 0.4
+
+            initial_condition = InitialCondition(; coordinates, mass, density)
+
+            function movement_function(initial_pos, t)
+                # Offset clamped particles to create stress
+                return initial_pos + SVector(0.0, t + 0.5)
+            end
+            is_moving(t) = true
+            prescribed_motion = PrescribedMotion(movement_function, is_moving)
+
+            # Create TLSPH system with energy calculator support
+            system = TotalLagrangianSPHSystem(initial_condition, smoothing_kernel,
+                                              smoothing_length, young_modulus,
+                                              poisson_ratio,
+                                              clamped_particles_motion=prescribed_motion,
+                                              n_clamped_particles=n_clamped_particles[i],
+                                              acceleration=(0.0, -2.0),
+                                              use_with_energy_calculator_callback=true)
+
+            semi = Semidiscretization(system)
+            ode = semidiscretize(semi, (0.0, 1.0))
+
+            # Create dummy ODE state vectors
+            v = zeros(2, TrixiParticles.n_integrated_particles(system))
+            u = coordinates[:, 1:TrixiParticles.n_integrated_particles(system)]
+            v_ode = vec(v)
+            u_ode = vec(u)
+
+            # Update system
+            TrixiParticles.initialize!(system, semi)
+            TrixiParticles.update_positions!(system, v, u, v_ode, u_ode, semi, 0.0)
+            TrixiParticles.update_quantities!(system, v, u, v_ode, u_ode, semi, 0.0)
+
+            # Set up test parameters
+            energy1 = Ref(0.0)
+            dt1 = 0.1
+
+            # Test that energy is integrated, i.e., values of instantaneous power
+            # are accumulated over time. This initial energy should just be an offset.
+            # Also, half the step size means half the energy increase.
+            energy2 = Ref(1.0)
+            dt2 = 0.05
+
+            TrixiParticles.update_energy_calculator!(energy1, v_ode, u_ode, system, semi,
+                                                     0.0, dt1)
+            TrixiParticles.update_energy_calculator!(energy2, v_ode, u_ode, system, semi,
+                                                     0.0, dt2)
+
+            if i == 1
+                @test isapprox(energy1[], 0.8)
+                @test isapprox(energy2[], 1.0 + 0.4)
+            elseif i == 2
+                # For very soft material, we can just pull up the top row of particles
+                # and the energy required is almost just the potential energy difference.
+                @test isapprox(energy1[], 0.4080357142857143)
+                @test isapprox(energy2[], 1.0 + 0.5 * 0.4080357142857143)
+            elseif i == 3
+                # For a stiffer material, the stress from the offset creates larger forces
+                # pulling the clamped particles back down, so we need a lot of energy
+                # to pull the material apart.
+                @test isapprox(energy1[], 803.9714285714281)
+                @test isapprox(energy2[], 1.0 + 0.5 * 803.9714285714281)
+            elseif i == 4
+                # For a very stiff material, the energy is even larger.
+                @test isapprox(energy1[], 80357.5428571428)
+                @test isapprox(energy2[], 1.0 + 0.5 * 80357.5428571428)
+            end
+        end
+    end
+end
diff --git a/test/examples/examples.jl b/test/examples/examples.jl
index 91ec13ef13..7a4bff9a38 100644
--- a/test/examples/examples.jl
+++ b/test/examples/examples.jl
@@ -41,6 +41,42 @@
             @test sol.retcode == ReturnCode.Success
             @test count_rhs_allocations(sol, semi) == 0
         end
+
+        @trixi_testset "structure/oscillating_beam_2d.jl with EnergyCalculatorCallback" begin
+            # Load variables
+            trixi_include(@__MODULE__,
+                          joinpath(examples_dir(), "structure", "oscillating_beam_2d.jl"),
+                          ode=nothing, sol=nothing)
+
+            # We simply clamp all particles, move them up against gravity, and verify that
+            # the energy calculated is just the potential energy difference.
+            movement_function(x, t) = x + SVector(0.0, t)
+            is_moving(t) = true
+            prescribed_motion = PrescribedMotion(movement_function, is_moving)
+
+            structure_system = TotalLagrangianSPHSystem(structure, smoothing_kernel,
+                                                        smoothing_length,
+                                                        material.E, material.nu,
+                                                        n_clamped_particles=nparticles(structure),
+                                                        acceleration=(0.0, -gravity),
+                                                        clamped_particles_motion=prescribed_motion,
+                                                        use_with_energy_calculator_callback=true)
+
+            semi = Semidiscretization(structure_system)
+            ode = semidiscretize(semi, (0.0, 1.0))
+
+            energy_calculator = EnergyCalculatorCallback{Float64}(; interval=1)
+
+            @trixi_test_nowarn sol = solve(ode, RDPK3SpFSAL49(), save_everystep=false,
+                                           callback=energy_calculator)
+
+            @test sol.retcode == ReturnCode.Success
+            @test count_rhs_allocations(sol, semi) == 0
+
+            # Potential energy difference should be m * g * h
+            @test isapprox(energy_calculator.affect!.energy[],
+                           sum(structure_system.mass) * gravity * 1)
+        end
     end
 
     @testset verbose=true "FSI" begin

From 9ebead743976b7b964fc15232997d5bf58ae007a Mon Sep 17 00:00:00 2001
From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com>
Date: Mon, 6 Oct 2025 15:00:15 +0200
Subject: [PATCH 06/16] Add function to retrieve calculated energy

---
 src/TrixiParticles.jl               |  2 +-
 src/callbacks/energy_calculator.jl  | 42 +++++++++++++++++++++++------
 test/callbacks/energy_calculator.jl |  1 +
 test/examples/examples.jl           |  2 +-
 4 files changed, 37 insertions(+), 10 deletions(-)

diff --git a/src/TrixiParticles.jl b/src/TrixiParticles.jl
index 70b733724a..3ba23449f4 100644
--- a/src/TrixiParticles.jl
+++ b/src/TrixiParticles.jl
@@ -70,7 +70,7 @@ export WeaklyCompressibleSPHSystem, EntropicallyDampedSPHSystem, TotalLagrangian
 export BoundaryZone, InFlow, OutFlow, BidirectionalFlow
 export InfoCallback, SolutionSavingCallback, DensityReinitializationCallback,
        PostprocessCallback, StepsizeCallback, UpdateCallback, SteadyStateReachedCallback,
-       SplitIntegrationCallback, EnergyCalculatorCallback
+       SplitIntegrationCallback, EnergyCalculatorCallback, calculated_energy
 export ContinuityDensity, SummationDensity
 export PenaltyForceGanzenmueller, TransportVelocityAdami, ParticleShiftingTechnique,
        ParticleShiftingTechniqueSun2017, ConsistentShiftingSun2019,
diff --git a/src/callbacks/energy_calculator.jl b/src/callbacks/energy_calculator.jl
index 7916b643c7..43bc882285 100644
--- a/src/callbacks/energy_calculator.jl
+++ b/src/callbacks/energy_calculator.jl
@@ -1,4 +1,3 @@
-
 """
     EnergyCalculatorCallback{ELTYPE}(; interval=1)
 
@@ -11,6 +10,8 @@ The energy is calculated by integrating the instantaneous power, which is the do
 of the force exerted by the clamped particle and its prescribed velocity.
 The callback uses a simple explicit Euler time integration for the energy calculation.
 
+The accumulated energy can be retrieved using [`calculated_energy`](@ref).
+
 !!! info "TLSPH System Configuration"
     The [`TotalLagrangianSPHSystem`](@ref) must be created with
     `use_with_energy_calculator_callback=true` to prepare the system for energy
@@ -42,17 +43,16 @@ system = TotalLagrangianSPHSystem(structure, smoothing_kernel, smoothing_length,
 # Create an energy calculator that runs every 2 time steps
 energy_cb = EnergyCalculatorCallback{eltype(system)}(; interval=2)
 
+# After the simulation, retrieve the calculated energy
+total_energy = calculated_energy(energy_cb)
+
 # output
-┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
-│ EnergyCalculatorCallback{Float64}                                                                │
-│ ═════════════════════════════════                                                                │
-│ interval: ……………………………………………………… 2                                                                │
-└──────────────────────────────────────────────────────────────────────────────────────────────────┘
+0.0
 ```
 """
 struct EnergyCalculatorCallback{T}
     interval :: Int
-    t        :: T
+    t        :: T # Time of last call
     energy   :: T
 end
 
@@ -63,6 +63,30 @@ function EnergyCalculatorCallback{ELTYPE}(; interval=1) where {ELTYPE <: Real}
     return DiscreteCallback(cb, cb, save_positions=(false, false))
 end
 
+"""
+    calculated_energy(cb::DiscreteCallback{<:Any, <:EnergyCalculatorCallback})
+
+Get the current calculated energy from an [`EnergyCalculatorCallback`](@ref).
+
+# Arguments
+- `cb`: The `DiscreteCallback` returned by [`EnergyCalculatorCallback`](@ref).
+
+# Examples
+```jldoctest; output = false
+# Before the simulation: create the callback
+energy_cb = EnergyCalculatorCallback{Float64}()
+
+# After the simulation: get the calculated energy
+total_energy = calculated_energy(energy_cb)
+
+# output
+0.0
+```
+"""
+function calculated_energy(cb::DiscreteCallback{<:Any, <:EnergyCalculatorCallback})
+    return cb.affect!.energy[]
+end
+
 # `condition`
 function (callback::EnergyCalculatorCallback)(u, t, integrator)
     (; interval) = callback
@@ -83,7 +107,9 @@ function (callback::EnergyCalculatorCallback)(integrator)
     v_ode, u_ode = integrator.u.x
     energy = callback.energy
 
-    update_energy_calculator!(energy, v_ode, u_ode, semi.systems[1], semi, t, dt)
+    foreach_system(semi) do system
+        update_energy_calculator!(energy, v_ode, u_ode, system, semi, t, dt)
+    end
 
     # Tell OrdinaryDiffEq that `u` has not been modified
     u_modified!(integrator, false)
diff --git a/test/callbacks/energy_calculator.jl b/test/callbacks/energy_calculator.jl
index 214efd9bbd..232f763b9f 100644
--- a/test/callbacks/energy_calculator.jl
+++ b/test/callbacks/energy_calculator.jl
@@ -5,6 +5,7 @@
         @test callback.affect!.interval == 1
         @test callback.affect!.t[] == 0.0
         @test callback.affect!.energy[] == 0.0
+        @test calculated_energy(callback) == 0.0
 
         # Test constructor with interval
         callback = EnergyCalculatorCallback{Float64}(; interval=5)
diff --git a/test/examples/examples.jl b/test/examples/examples.jl
index 7a4bff9a38..3222f044a3 100644
--- a/test/examples/examples.jl
+++ b/test/examples/examples.jl
@@ -74,7 +74,7 @@
             @test count_rhs_allocations(sol, semi) == 0
 
             # Potential energy difference should be m * g * h
-            @test isapprox(energy_calculator.affect!.energy[],
+            @test isapprox(calculated_energy(energy_calculator),
                            sum(structure_system.mass) * gravity * 1)
         end
     end

From 7298243e03d8bfb316f4af67012856d052ed7c32 Mon Sep 17 00:00:00 2001
From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com>
Date: Mon, 6 Oct 2025 16:17:39 +0200
Subject: [PATCH 07/16] Fix for split integration

---
 src/callbacks/energy_calculator.jl | 5 +++++
 1 file changed, 5 insertions(+)

diff --git a/src/callbacks/energy_calculator.jl b/src/callbacks/energy_calculator.jl
index 43bc882285..ca911cc10c 100644
--- a/src/callbacks/energy_calculator.jl
+++ b/src/callbacks/energy_calculator.jl
@@ -117,6 +117,10 @@ function (callback::EnergyCalculatorCallback)(integrator)
     return integrator
 end
 
+function update_energy_calculator!(energy, v_ode, u_ode, system, semi, t, dt)
+    return energy
+end
+
 function update_energy_calculator!(energy, v_ode, u_ode,
                                    system::AbstractStructureSystem, semi, t, dt)
     @trixi_timeit timer() "calculate energy" begin
@@ -142,6 +146,7 @@ function update_energy_calculator!(energy, v_ode, u_ode,
 
             interact!(dv_combined, v, u, v_neighbor, u_neighbor,
                       system, neighbor_system, semi,
+                      integrate_tlsph=true, # Required when using split integration
                       eachparticle=eachparticle)
         end
 

From b8594d551b16c3a24bd7c11853090a3e814441d4 Mon Sep 17 00:00:00 2001
From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com>
Date: Mon, 6 Oct 2025 16:28:02 +0200
Subject: [PATCH 08/16] Use `OffsetMatrix` instead of `CombinedMatrix` to avoid
 allocations

---
 src/callbacks/energy_calculator.jl | 72 ++++++++++++++----------------
 1 file changed, 33 insertions(+), 39 deletions(-)

diff --git a/src/callbacks/energy_calculator.jl b/src/callbacks/energy_calculator.jl
index ca911cc10c..290be9c298 100644
--- a/src/callbacks/energy_calculator.jl
+++ b/src/callbacks/energy_calculator.jl
@@ -134,8 +134,12 @@ function update_energy_calculator!(energy, v_ode, u_ode,
         dv_clamped = system.cache.dv_clamped
         set_zero!(dv_clamped)
 
-        dv = zero(wrap_v(v_ode, system, semi))
-        dv_combined = CombinedMatrix(dv, dv_clamped)
+        # Create a view of `dv_clamped` that can be indexed as
+        # `(n_integrated_particles(system) + 1):nparticles(system)`.
+        # We pass this to `interact!` below so that it can compute
+        # the forces on the clamped particles without having to allocate
+        # a large matrix for all particles.
+        dv_offset = OffsetMatrix(dv_clamped, n_integrated_particles(system))
         v = wrap_v(v_ode, system, semi)
         u = wrap_u(u_ode, system, semi)
         eachparticle = (n_integrated_particles(system) + 1):nparticles(system)
@@ -144,23 +148,21 @@ function update_energy_calculator!(energy, v_ode, u_ode,
             v_neighbor = wrap_v(v_ode, neighbor_system, semi)
             u_neighbor = wrap_u(u_ode, neighbor_system, semi)
 
-            interact!(dv_combined, v, u, v_neighbor, u_neighbor,
+            interact!(dv_offset, v, u, v_neighbor, u_neighbor,
                       system, neighbor_system, semi,
                       integrate_tlsph=true, # Required when using split integration
                       eachparticle=eachparticle)
         end
 
         @threaded semi for particle in eachparticle
-            add_acceleration!(dv_combined, particle, system)
-            add_source_terms_inner!(dv_combined, v, u, particle, system,
+            add_acceleration!(dv_offset, particle, system)
+            add_source_terms_inner!(dv_offset, v, u, particle, system,
                                     source_terms(system), t)
         end
 
-        n_clamped_particles = nparticles(system) - n_integrated_particles(system)
-        for clamped_particle in 1:n_clamped_particles
-            particle = clamped_particle + n_integrated_particles(system)
+        for particle in (n_integrated_particles(system) + 1):nparticles(system)
             velocity = current_velocity(nothing, system, particle)
-            dv_particle = extract_svector(dv_clamped, system, clamped_particle)
+            dv_particle = extract_svector(dv_offset, system, particle)
 
             # The force on the clamped particle is mass times acceleration
             F_particle = system.mass[particle] * dv_particle
@@ -196,43 +198,35 @@ function Base.show(io::IO, ::MIME"text/plain",
     end
 end
 
-# Data type that combines two matrices to behave like a single matrix.
-# This is used above to combine the `dv` for integrated particles and the `dv_clamped`
-# for clamped particles into a single matrix that can be passed to `interact!`.
-struct CombinedMatrix{T, M1, M2} <: AbstractMatrix{T}
-    matrix1::M1
-    matrix2::M2
-
-    function CombinedMatrix(matrix1, matrix2)
-        @assert size(matrix1, 1) == size(matrix2, 1)
-        @assert eltype(matrix1) == eltype(matrix2)
-
-        new{eltype(matrix1), typeof(matrix1), typeof(matrix2)}(matrix1, matrix2)
+# Data type that represents a matrix with an offset in the second dimension.
+# For example, `om = OffsetMatrix(A, offset)` starts at `om[:, offset + 1]`, which is
+# the same as `A[:, 1]` and it ends at `om[:, offset + size(A, 2)]`, which is
+# the same as `A[:, end]`.
+# This is used above to compute the acceleration only for clamped particles (which are the
+# last particles in the system) without having to pass a large matrix for all particles.
+struct OffsetMatrix{T, M} <: AbstractMatrix{T}
+    matrix::M
+    offset::Int
+
+    function OffsetMatrix(matrix, offset)
+        new{eltype(matrix), typeof(matrix)}(matrix, offset)
     end
 end
 
-@inline function Base.size(cm::CombinedMatrix)
-    return (size(cm.matrix1, 1), size(cm.matrix1, 2) + size(cm.matrix2, 2))
+@inline function Base.size(om::OffsetMatrix)
+    return (size(om.matrix, 1), size(om.matrix, 2) + om.offset)
 end
 
-@inline function Base.getindex(cm::CombinedMatrix, i, j)
-    @boundscheck checkbounds(cm, i, j)
+@inline function Base.getindex(om::OffsetMatrix, i, j)
+    @boundscheck checkbounds(om, i, j)
+    @boundscheck j > om.offset
 
-    length1 = size(cm.matrix1, 2)
-    if j <= length1
-        return @inbounds cm.matrix1[i, j]
-    else
-        return @inbounds cm.matrix2[i, j - length1]
-    end
+    return @inbounds om.matrix[i, j - om.offset]
 end
 
-@inline function Base.setindex!(cm::CombinedMatrix, value, i, j)
-    @boundscheck checkbounds(cm, i, j)
+@inline function Base.setindex!(om::OffsetMatrix, value, i, j)
+    @boundscheck checkbounds(om, i, j)
+    @boundscheck j > om.offset
 
-    length1 = size(cm.matrix1, 2)
-    if j <= length1
-        return @inbounds cm.matrix1[i, j] = value
-    else
-        return @inbounds cm.matrix2[i, j - length1] = value
-    end
+    return @inbounds om.matrix[i, j - om.offset] = value
 end

From 6016e450034cc999dce92708b29b101e1d0a867a Mon Sep 17 00:00:00 2001
From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com>
Date: Mon, 6 Oct 2025 16:28:33 +0200
Subject: [PATCH 09/16] Fix return value

---
 src/callbacks/energy_calculator.jl | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/src/callbacks/energy_calculator.jl b/src/callbacks/energy_calculator.jl
index 290be9c298..ed996633ac 100644
--- a/src/callbacks/energy_calculator.jl
+++ b/src/callbacks/energy_calculator.jl
@@ -173,6 +173,8 @@ function update_energy_calculator!(energy, v_ode, u_ode,
             energy[] -= dot(F_particle, velocity) * dt
         end
     end
+
+    return energy
 end
 
 function Base.show(io::IO, cb::DiscreteCallback{<:Any, <:EnergyCalculatorCallback})

From 93273d92b1d075428ae72aa20db57f4e5b951d1c Mon Sep 17 00:00:00 2001
From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com>
Date: Mon, 6 Oct 2025 16:36:30 +0200
Subject: [PATCH 10/16] Fix tests

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

diff --git a/test/examples/examples.jl b/test/examples/examples.jl
index 3222f044a3..589fa2dfa5 100644
--- a/test/examples/examples.jl
+++ b/test/examples/examples.jl
@@ -67,7 +67,7 @@
 
             energy_calculator = EnergyCalculatorCallback{Float64}(; interval=1)
 
-            @trixi_test_nowarn sol = solve(ode, RDPK3SpFSAL49(), save_everystep=false,
+            sol = @trixi_test_nowarn solve(ode, RDPK3SpFSAL49(), save_everystep=false,
                                            callback=energy_calculator)
 
             @test sol.retcode == ReturnCode.Success

From 6ae7dc1340c808f4598bcf2158cd518209834282 Mon Sep 17 00:00:00 2001
From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com>
Date: Wed, 15 Oct 2025 12:16:21 +0200
Subject: [PATCH 11/16] Only calculate energy for a single system

---
 src/callbacks/energy_calculator.jl  | 100 +++++++++-------------------
 src/general/semidiscretization.jl   |   2 +-
 test/callbacks/energy_calculator.jl |  35 +++++++---
 3 files changed, 60 insertions(+), 77 deletions(-)

diff --git a/src/callbacks/energy_calculator.jl b/src/callbacks/energy_calculator.jl
index ed996633ac..f71620be06 100644
--- a/src/callbacks/energy_calculator.jl
+++ b/src/callbacks/energy_calculator.jl
@@ -50,14 +50,24 @@ total_energy = calculated_energy(energy_cb)
 0.0
 ```
 """
-struct EnergyCalculatorCallback{T}
-    interval :: Int
-    t        :: T # Time of last call
-    energy   :: T
+struct EnergyCalculatorCallback{T, DV, EP}
+    interval     :: Int
+    t            :: T # Time of last call
+    energy       :: T
+    system_index :: Int
+    dv           :: DV
+    eachparticle :: EP
 end
 
-function EnergyCalculatorCallback{ELTYPE}(; interval=1) where {ELTYPE <: Real}
-    cb = EnergyCalculatorCallback(interval, Ref(zero(ELTYPE)), Ref(zero(ELTYPE)))
+function EnergyCalculatorCallback{ELTYPE}(system, semi; interval=1,
+                                          eachparticle=eachparticle(system)) where {ELTYPE <: Real}
+    system_index = system_indices(system, semi)
+
+    # Allocate buffer to write accelerations for all particles (including clamped ones)
+    dv = allocate(semi.parallelization_backend, ELTYPE, (ndims(system), nparticles(system)))
+
+    cb = EnergyCalculatorCallback(interval, Ref(zero(ELTYPE)), Ref(zero(ELTYPE)),
+                                  system_index, dv, eachparticle)
 
     # The first one is the `condition`, the second the `affect!`
     return DiscreteCallback(cb, cb, save_positions=(false, false))
@@ -107,9 +117,9 @@ function (callback::EnergyCalculatorCallback)(integrator)
     v_ode, u_ode = integrator.u.x
     energy = callback.energy
 
-    foreach_system(semi) do system
-        update_energy_calculator!(energy, v_ode, u_ode, system, semi, t, dt)
-    end
+    system = semi.systems[callback.system_index]
+    update_energy_calculator!(energy, system, callback.eachparticle, callback.dv,
+                              v_ode, u_ode, semi, t, dt)
 
     # Tell OrdinaryDiffEq that `u` has not been modified
     u_modified!(integrator, false)
@@ -117,12 +127,8 @@ function (callback::EnergyCalculatorCallback)(integrator)
     return integrator
 end
 
-function update_energy_calculator!(energy, v_ode, u_ode, system, semi, t, dt)
-    return energy
-end
-
-function update_energy_calculator!(energy, v_ode, u_ode,
-                                   system::AbstractStructureSystem, semi, t, dt)
+function update_energy_calculator!(energy, system, eachparticle, dv,
+                                   v_ode, u_ode, semi, t, dt)
     @trixi_timeit timer() "calculate energy" 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.
@@ -131,46 +137,37 @@ function update_energy_calculator!(energy, v_ode, u_ode,
             @notimeit timer() update_systems_and_nhs(v_ode, u_ode, semi, t)
         end
 
-        dv_clamped = system.cache.dv_clamped
-        set_zero!(dv_clamped)
+        set_zero!(dv)
 
-        # Create a view of `dv_clamped` that can be indexed as
-        # `(n_integrated_particles(system) + 1):nparticles(system)`.
-        # We pass this to `interact!` below so that it can compute
-        # the forces on the clamped particles without having to allocate
-        # a large matrix for all particles.
-        dv_offset = OffsetMatrix(dv_clamped, n_integrated_particles(system))
         v = wrap_v(v_ode, system, semi)
         u = wrap_u(u_ode, system, semi)
-        eachparticle = (n_integrated_particles(system) + 1):nparticles(system)
 
         foreach_system(semi) do neighbor_system
             v_neighbor = wrap_v(v_ode, neighbor_system, semi)
             u_neighbor = wrap_u(u_ode, neighbor_system, semi)
 
-            interact!(dv_offset, v, u, v_neighbor, u_neighbor,
+            interact!(dv, v, u, v_neighbor, u_neighbor,
                       system, neighbor_system, semi,
                       integrate_tlsph=true, # Required when using split integration
                       eachparticle=eachparticle)
         end
 
         @threaded semi for particle in eachparticle
-            add_acceleration!(dv_offset, particle, system)
-            add_source_terms_inner!(dv_offset, v, u, particle, system,
-                                    source_terms(system), t)
+            add_acceleration!(dv, particle, system)
+            add_source_terms_inner!(dv, v, u, particle, system, source_terms(system), t)
         end
 
-        for particle in (n_integrated_particles(system) + 1):nparticles(system)
-            velocity = current_velocity(nothing, system, particle)
-            dv_particle = extract_svector(dv_offset, system, particle)
+        for particle in eachparticle
+            velocity = current_velocity(v, system, particle)
+            dv_particle = extract_svector(dv, system, particle)
 
             # The force on the clamped particle is mass times acceleration
             F_particle = system.mass[particle] * dv_particle
 
             # To obtain energy, we need to integrate the instantaneous power.
-            # Instantaneous power is force done BY the particle times prescribed velocity.
-            # The work done BY the particle is the negative of the work done ON it.
-            energy[] -= dot(F_particle, velocity) * dt
+            # Instantaneous power is force applied BY the particle times its velocity.
+            # The force applied BY the particle is the negative of the force applied ON it.
+            energy[] += dot(-F_particle, velocity) * dt
         end
     end
 
@@ -185,7 +182,7 @@ function Base.show(io::IO, cb::DiscreteCallback{<:Any, <:EnergyCalculatorCallbac
 end
 
 function Base.show(io::IO, ::MIME"text/plain",
-                   cb::DiscreteCallback{<:Any, EnergyCalculatorCallback{T}}) where {T}
+                   cb::DiscreteCallback{<:Any, <:EnergyCalculatorCallback})
     @nospecialize cb # reduce precompilation time
 
     if get(io, :compact, false)
@@ -199,36 +196,3 @@ function Base.show(io::IO, ::MIME"text/plain",
         summary_box(io, "EnergyCalculatorCallback{$ELTYPE}", setup)
     end
 end
-
-# Data type that represents a matrix with an offset in the second dimension.
-# For example, `om = OffsetMatrix(A, offset)` starts at `om[:, offset + 1]`, which is
-# the same as `A[:, 1]` and it ends at `om[:, offset + size(A, 2)]`, which is
-# the same as `A[:, end]`.
-# This is used above to compute the acceleration only for clamped particles (which are the
-# last particles in the system) without having to pass a large matrix for all particles.
-struct OffsetMatrix{T, M} <: AbstractMatrix{T}
-    matrix::M
-    offset::Int
-
-    function OffsetMatrix(matrix, offset)
-        new{eltype(matrix), typeof(matrix)}(matrix, offset)
-    end
-end
-
-@inline function Base.size(om::OffsetMatrix)
-    return (size(om.matrix, 1), size(om.matrix, 2) + om.offset)
-end
-
-@inline function Base.getindex(om::OffsetMatrix, i, j)
-    @boundscheck checkbounds(om, i, j)
-    @boundscheck j > om.offset
-
-    return @inbounds om.matrix[i, j - om.offset]
-end
-
-@inline function Base.setindex!(om::OffsetMatrix, value, i, j)
-    @boundscheck checkbounds(om, i, j)
-    @boundscheck j > om.offset
-
-    return @inbounds om.matrix[i, j - om.offset] = value
-end
diff --git a/src/general/semidiscretization.jl b/src/general/semidiscretization.jl
index ceba8b5bdf..4c3206528f 100644
--- a/src/general/semidiscretization.jl
+++ b/src/general/semidiscretization.jl
@@ -232,7 +232,7 @@ end
 @inline function system_indices(system, semi)
     # Note that this takes only about 5 ns, while mapping systems to indices with a `Dict`
     # is ~30x slower because `hash(::System)` is very slow.
-    index = findfirst(==(system), semi.systems)
+    index = findfirst(s -> s === system, semi.systems)
 
     if isnothing(index)
         throw(ArgumentError("system is not in the semidiscretization"))
diff --git a/test/callbacks/energy_calculator.jl b/test/callbacks/energy_calculator.jl
index 232f763b9f..06a6b5285a 100644
--- a/test/callbacks/energy_calculator.jl
+++ b/test/callbacks/energy_calculator.jl
@@ -1,24 +1,40 @@
 @testset verbose=true "EnergyCalculatorCallback" begin
+    # Mock system
+    struct MockSystem <: TrixiParticles.AbstractSystem{2} end
+    TrixiParticles.nparticles(::MockSystem) = 4
+
+    function TrixiParticles.create_neighborhood_search(neighborhood_search,
+                                                       system::MockSystem, neighbor)
+        return nothing
+    end
+
+    system = MockSystem()
+    semi = Semidiscretization(system)
+
     @testset "Constructor and Basic Properties" begin
         # Test default constructor
-        callback = EnergyCalculatorCallback{Float64}()
+        callback = EnergyCalculatorCallback{Float64}(system, semi)
+        @test callback.affect!.system_index == 1
         @test callback.affect!.interval == 1
         @test callback.affect!.t[] == 0.0
         @test callback.affect!.energy[] == 0.0
+        @test callback.affect!.dv isa Array{Float64, 2}
+        @test size(callback.affect!.dv) == (2, 4)
+        @test callback.affect!.eachparticle == 1:4
         @test calculated_energy(callback) == 0.0
 
         # Test constructor with interval
-        callback = EnergyCalculatorCallback{Float64}(; interval=5)
+        callback = EnergyCalculatorCallback{Float64}(system, semi; interval=5)
         @test callback.affect!.interval == 5
 
         # Test with specific element type
-        callback = EnergyCalculatorCallback{Float32}(; interval=2)
+        callback = EnergyCalculatorCallback{Float32}(system, semi; interval=2)
         @test eltype(callback.affect!.energy) == Float32
         @test eltype(callback.affect!.t) == Float32
     end
 
     @testset "show" begin
-        callback = EnergyCalculatorCallback{Float64}(; interval=10)
+        callback = EnergyCalculatorCallback{Float64}(system, semi; interval=10)
 
         # Test compact representation
         show_compact = "EnergyCalculatorCallback{Float64}(interval=10)"
@@ -96,10 +112,13 @@
             energy2 = Ref(1.0)
             dt2 = 0.05
 
-            TrixiParticles.update_energy_calculator!(energy1, v_ode, u_ode, system, semi,
-                                                     0.0, dt1)
-            TrixiParticles.update_energy_calculator!(energy2, v_ode, u_ode, system, semi,
-                                                     0.0, dt2)
+            eachparticle = (TrixiParticles.n_integrated_particles(system) + 1):nparticles(system)
+            dv = zeros(2, nparticles(system))
+
+            TrixiParticles.update_energy_calculator!(energy1, system, eachparticle, dv,
+                                                     v_ode, u_ode, semi, 0.0, dt1)
+            TrixiParticles.update_energy_calculator!(energy2, system, eachparticle, dv,
+                                                     v_ode, u_ode, semi, 0.0, dt2)
 
             if i == 1
                 @test isapprox(energy1[], 0.8)

From 8fd43f9462f1916e2a5e475ce6f8ee3e96533bab Mon Sep 17 00:00:00 2001
From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com>
Date: Wed, 15 Oct 2025 12:27:53 +0200
Subject: [PATCH 12/16] Implement `only_compute_force_on_fluid`

---
 src/callbacks/energy_calculator.jl  | 38 +++++++++++++++++++----------
 test/callbacks/energy_calculator.jl | 16 +++++++++---
 2 files changed, 37 insertions(+), 17 deletions(-)

diff --git a/src/callbacks/energy_calculator.jl b/src/callbacks/energy_calculator.jl
index f71620be06..55ba8114af 100644
--- a/src/callbacks/energy_calculator.jl
+++ b/src/callbacks/energy_calculator.jl
@@ -51,23 +51,26 @@ total_energy = calculated_energy(energy_cb)
 ```
 """
 struct EnergyCalculatorCallback{T, DV, EP}
-    interval     :: Int
-    t            :: T # Time of last call
-    energy       :: T
-    system_index :: Int
-    dv           :: DV
-    eachparticle :: EP
+    interval                      :: Int
+    t                             :: T # Time of last call
+    energy                        :: T
+    system_index                  :: Int
+    dv                            :: DV
+    eachparticle                  :: EP
+    only_compute_force_on_fluid   :: Bool
 end
 
 function EnergyCalculatorCallback{ELTYPE}(system, semi; interval=1,
-                                          eachparticle=eachparticle(system)) where {ELTYPE <: Real}
+                                          eachparticle=eachparticle(system),
+                                          only_compute_force_on_fluid=false) where {ELTYPE <: Real}
     system_index = system_indices(system, semi)
 
     # Allocate buffer to write accelerations for all particles (including clamped ones)
     dv = allocate(semi.parallelization_backend, ELTYPE, (ndims(system), nparticles(system)))
 
     cb = EnergyCalculatorCallback(interval, Ref(zero(ELTYPE)), Ref(zero(ELTYPE)),
-                                  system_index, dv, eachparticle)
+                                  system_index, dv, eachparticle,
+                                  only_compute_force_on_fluid)
 
     # The first one is the `condition`, the second the `affect!`
     return DiscreteCallback(cb, cb, save_positions=(false, false))
@@ -118,7 +121,8 @@ function (callback::EnergyCalculatorCallback)(integrator)
     energy = callback.energy
 
     system = semi.systems[callback.system_index]
-    update_energy_calculator!(energy, system, callback.eachparticle, callback.dv,
+    update_energy_calculator!(energy, system, callback.eachparticle,
+                              callback.only_compute_force_on_fluid, callback.dv,
                               v_ode, u_ode, semi, t, dt)
 
     # Tell OrdinaryDiffEq that `u` has not been modified
@@ -127,7 +131,8 @@ function (callback::EnergyCalculatorCallback)(integrator)
     return integrator
 end
 
-function update_energy_calculator!(energy, system, eachparticle, dv,
+function update_energy_calculator!(energy, system, eachparticle,
+                                   only_compute_force_on_fluid, dv,
                                    v_ode, u_ode, semi, t, dt)
     @trixi_timeit timer() "calculate energy" begin
         # Update quantities that are stored in the systems. These quantities (e.g. pressure)
@@ -143,6 +148,11 @@ function update_energy_calculator!(energy, system, eachparticle, dv,
         u = wrap_u(u_ode, system, semi)
 
         foreach_system(semi) do neighbor_system
+            if only_compute_force_on_fluid && !(neighbor_system isa AbstractFluidSystem)
+                # Not a fluid system, ignore this system
+                return
+            end
+
             v_neighbor = wrap_v(v_ode, neighbor_system, semi)
             u_neighbor = wrap_u(u_ode, neighbor_system, semi)
 
@@ -152,9 +162,11 @@ function update_energy_calculator!(energy, system, eachparticle, dv,
                       eachparticle=eachparticle)
         end
 
-        @threaded semi for particle in eachparticle
-            add_acceleration!(dv, particle, system)
-            add_source_terms_inner!(dv, v, u, particle, system, source_terms(system), t)
+        if !only_compute_force_on_fluid
+            @threaded semi for particle in eachparticle
+                add_acceleration!(dv, particle, system)
+                add_source_terms_inner!(dv, v, u, particle, system, source_terms(system), t)
+            end
         end
 
         for particle in eachparticle
diff --git a/test/callbacks/energy_calculator.jl b/test/callbacks/energy_calculator.jl
index 06a6b5285a..74be5dc1b5 100644
--- a/test/callbacks/energy_calculator.jl
+++ b/test/callbacks/energy_calculator.jl
@@ -112,13 +112,19 @@
             energy2 = Ref(1.0)
             dt2 = 0.05
 
+            # Test `only_compute_force_on_fluid`
+            energy3 = Ref(0.0)
+            dt3 = 0.1
+
             eachparticle = (TrixiParticles.n_integrated_particles(system) + 1):nparticles(system)
             dv = zeros(2, nparticles(system))
 
-            TrixiParticles.update_energy_calculator!(energy1, system, eachparticle, dv,
-                                                     v_ode, u_ode, semi, 0.0, dt1)
-            TrixiParticles.update_energy_calculator!(energy2, system, eachparticle, dv,
-                                                     v_ode, u_ode, semi, 0.0, dt2)
+            TrixiParticles.update_energy_calculator!(energy1, system, eachparticle, false,
+                                                     dv, v_ode, u_ode, semi, 0.0, dt1)
+            TrixiParticles.update_energy_calculator!(energy2, system, eachparticle, false,
+                                                     dv, v_ode, u_ode, semi, 0.0, dt2)
+            TrixiParticles.update_energy_calculator!(energy3, system, eachparticle, true,
+                                                     dv, v_ode, u_ode, semi, 0.0, dt3)
 
             if i == 1
                 @test isapprox(energy1[], 0.8)
@@ -139,6 +145,8 @@
                 @test isapprox(energy1[], 80357.5428571428)
                 @test isapprox(energy2[], 1.0 + 0.5 * 80357.5428571428)
             end
+
+            @test isapprox(energy3[], 0.0)
         end
     end
 end

From 4e87fca5708af4b6306381aa4826abdb24658ee0 Mon Sep 17 00:00:00 2001
From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com>
Date: Wed, 15 Oct 2025 12:36:14 +0200
Subject: [PATCH 13/16] Fix example test

---
 test/examples/examples.jl | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/test/examples/examples.jl b/test/examples/examples.jl
index 589fa2dfa5..192ef7cfdb 100644
--- a/test/examples/examples.jl
+++ b/test/examples/examples.jl
@@ -65,7 +65,8 @@
             semi = Semidiscretization(structure_system)
             ode = semidiscretize(semi, (0.0, 1.0))
 
-            energy_calculator = EnergyCalculatorCallback{Float64}(; interval=1)
+            energy_calculator = EnergyCalculatorCallback{Float64}(structure_system, semi;
+                                                                  interval=1)
 
             sol = @trixi_test_nowarn solve(ode, RDPK3SpFSAL49(), save_everystep=false,
                                            callback=energy_calculator)

From 019be6f54ecce89ded14f8bae6bc6c2d45cf8ee1 Mon Sep 17 00:00:00 2001
From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com>
Date: Wed, 15 Oct 2025 15:33:45 +0200
Subject: [PATCH 14/16] Add second test with FSI

---
 test/examples/examples.jl | 78 +++++++++++++++++++++++++++++++++++++--
 1 file changed, 75 insertions(+), 3 deletions(-)

diff --git a/test/examples/examples.jl b/test/examples/examples.jl
index 192ef7cfdb..8b17226d88 100644
--- a/test/examples/examples.jl
+++ b/test/examples/examples.jl
@@ -43,7 +43,7 @@
         end
 
         @trixi_testset "structure/oscillating_beam_2d.jl with EnergyCalculatorCallback" begin
-            # Load variables
+            # Load variables from the example
             trixi_include(@__MODULE__,
                           joinpath(examples_dir(), "structure", "oscillating_beam_2d.jl"),
                           ode=nothing, sol=nothing)
@@ -59,8 +59,7 @@
                                                         material.E, material.nu,
                                                         n_clamped_particles=nparticles(structure),
                                                         acceleration=(0.0, -gravity),
-                                                        clamped_particles_motion=prescribed_motion,
-                                                        use_with_energy_calculator_callback=true)
+                                                        clamped_particles_motion=prescribed_motion)
 
             semi = Semidiscretization(structure_system)
             ode = semidiscretize(semi, (0.0, 1.0))
@@ -81,6 +80,79 @@
     end
 
     @testset verbose=true "FSI" begin
+        @trixi_testset "fluid/hydrostatic_water_column_2d.jl with moving TLSPH walls" begin
+            # In this test, we move a water-filled tank up against gravity by 1 unit
+            # and verify that the energy calculated by the `EnergyCalculatorCallback`
+            # matches the expected potential energy.
+
+            # Load variables from the example
+            @trixi_test_nowarn trixi_include(@__MODULE__,
+                                             joinpath(examples_dir(), "fluid",
+                                                      "hydrostatic_water_column_2d.jl"),
+                                             # Use tank without airspace
+                                             tank_size=(1.0, 0.9),
+                                             # Higher speed of sound for smaller error
+                                             # due to compressibility.
+                                             sound_speed=50.0,
+                                             ode=nothing, sol=nothing)
+
+            # Move the tank up against gravity by 1 unit over 1 second
+            function movement_function(x, t)
+                return x + SVector(0.0, 0.5 * sin(pi * (t - 0.5)) + 0.5)
+            end
+            is_moving(t) = true
+            prescribed_motion = PrescribedMotion(movement_function, is_moving)
+
+            # Create TLSPH system for the tank walls and clamp all particles.
+            # This is identical to a `WallBoundarySystem`, but now we can
+            # use the `EnergyCalculatorCallback` to compute the energy.
+            boundary_spacing = tank.boundary.particle_spacing
+            tlsph_kernel = WendlandC2Kernel{2}()
+            tlsph_smoothing_length = sqrt(2) * boundary_spacing
+
+            tlsph_system = TotalLagrangianSPHSystem(tank.boundary, tlsph_kernel,
+                                                    tlsph_smoothing_length,
+                                                    1.0e6, 0.3,
+                                                    n_clamped_particles=nparticles(tank.boundary),
+                                                    acceleration=(0.0, -gravity),
+                                                    clamped_particles_motion=prescribed_motion,
+                                                    boundary_model=boundary_model)
+
+            semi = Semidiscretization(fluid_system, tlsph_system,
+                                      parallelization_backend=PolyesterBackend())
+            ode = semidiscretize(semi, (0.0, 1.0))
+
+            # Energy calculators for fluid + tank and fluid only
+            energy_calculator1 = EnergyCalculatorCallback{eltype(tlsph_system)}(tlsph_system,
+                                                                                semi;
+                                                                                interval=1)
+            energy_calculator2 = EnergyCalculatorCallback{eltype(tlsph_system)}(tlsph_system,
+                                                                                semi;
+                                                                                interval=1,
+                                                                                only_compute_force_on_fluid=true)
+
+            sol = @trixi_test_nowarn solve(ode, RDPK3SpFSAL35(), save_everystep=false,
+                                           callback=CallbackSet(info_callback,
+                                                                energy_calculator1,
+                                                                energy_calculator2))
+
+            @test sol.retcode == ReturnCode.Success
+            @test count_rhs_allocations(sol, semi) == 0
+
+            # Potential energy difference should be m * g * h and h = 1
+            expected_energy_fluid = sum(fluid_system.mass) * gravity * 1
+            expected_energy_tank = sum(tlsph_system.mass) * gravity * 1
+
+            # We don't expect very accurate results here because the fluid is weakly
+            # compressible and is deformed during the simulation.
+            # A slower prescribed motion (e.g., over 2 seconds instead of 1) or a higher
+            # speed of sound in the fluid would improve accuracy (and increase runtime).
+            @test isapprox(calculated_energy(energy_calculator1),
+                           expected_energy_fluid + expected_energy_tank, rtol=5e-4)
+            @test isapprox(calculated_energy(energy_calculator2), expected_energy_fluid,
+                           rtol=5e-4)
+        end
+
         @trixi_testset "fsi/falling_water_column_2d.jl" begin
             @trixi_test_nowarn trixi_include(@__MODULE__,
                                              joinpath(examples_dir(), "fsi",

From 96e3469a8c3f2503c90d58d9ad400220574ce4b9 Mon Sep 17 00:00:00 2001
From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com>
Date: Wed, 15 Oct 2025 15:35:34 +0200
Subject: [PATCH 15/16] Undo changes in TLSPH system

---
 .../structure/total_lagrangian_sph/system.jl       | 14 +-------------
 1 file changed, 1 insertion(+), 13 deletions(-)

diff --git a/src/schemes/structure/total_lagrangian_sph/system.jl b/src/schemes/structure/total_lagrangian_sph/system.jl
index c88bd42534..0900ac42b9 100644
--- a/src/schemes/structure/total_lagrangian_sph/system.jl
+++ b/src/schemes/structure/total_lagrangian_sph/system.jl
@@ -91,7 +91,6 @@ end
 function TotalLagrangianSPHSystem(initial_condition, smoothing_kernel, smoothing_length,
                                   young_modulus, poisson_ratio;
                                   n_clamped_particles=0, clamped_particles_motion=nothing,
-                                  use_with_energy_calculator_callback=false,
                                   boundary_model=nothing,
                                   acceleration=ntuple(_ -> 0.0,
                                                       ndims(smoothing_kernel)),
@@ -129,9 +128,7 @@ function TotalLagrangianSPHSystem(initial_condition, smoothing_kernel, smoothing
     initialize_prescribed_motion!(clamped_particles_motion, initial_condition,
                                   n_clamped_particles)
 
-    cache = (; create_cache_tlsph(clamped_particles_motion, initial_condition)...,
-             create_cache_tlsph(Val(use_with_energy_calculator_callback),
-                                initial_condition, n_clamped_particles)...)
+    cache = create_cache_tlsph(clamped_particles_motion, initial_condition)
 
     return TotalLagrangianSPHSystem(initial_condition, initial_coordinates,
                                     current_coordinates, mass, correction_matrix,
@@ -152,15 +149,6 @@ function create_cache_tlsph(::PrescribedMotion, initial_condition)
     return (; velocity, acceleration)
 end
 
-create_cache_tlsph(::Val{false}, initial_condition, n_clamped_particles) = (;)
-
-function create_cache_tlsph(::Val{true}, initial_condition, n_clamped_particles)
-    dv_clamped = Array{eltype(initial_condition)}(undef, ndims(initial_condition),
-                                                  n_clamped_particles)
-
-    return (; dv_clamped)
-end
-
 @inline function Base.eltype(::TotalLagrangianSPHSystem{<:Any, <:Any, ELTYPE}) where {ELTYPE}
     return ELTYPE
 end

From c9df08f0e4ada822a1350a991630e868c13b81af Mon Sep 17 00:00:00 2001
From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com>
Date: Wed, 15 Oct 2025 16:31:16 +0200
Subject: [PATCH 16/16] Update docs

---
 src/callbacks/energy_calculator.jl | 78 ++++++++++++++++--------------
 1 file changed, 42 insertions(+), 36 deletions(-)

diff --git a/src/callbacks/energy_calculator.jl b/src/callbacks/energy_calculator.jl
index 55ba8114af..9f5e9c917b 100644
--- a/src/callbacks/energy_calculator.jl
+++ b/src/callbacks/energy_calculator.jl
@@ -1,47 +1,52 @@
 """
-    EnergyCalculatorCallback{ELTYPE}(; interval=1)
+    EnergyCalculatorCallback(system, semi; interval=1,
+                             eachparticle=(n_integrated_particles(system) + 1):nparticles(system),
+                             only_compute_force_on_fluid=false)
 
-Callback to calculate the energy contribution from clamped particles in
-[`TotalLagrangianSPHSystem`](@ref)s.
+Callback to integrate the work done by particles in a
+[`TotalLagrangianSPHSystem`](@ref).
 
-This callback computes the work done by clamped particles that are moving according to a
-[`PrescribedMotion`](@ref).
-The energy is calculated by integrating the instantaneous power, which is the dot product
-of the force exerted by the clamped particle and its prescribed velocity.
-The callback uses a simple explicit Euler time integration for the energy calculation.
+With the default arguments it tracks the energy contribution of the clamped particles
+that follow a [`PrescribedMotion`](@ref). By selecting different particles it can also be
+used to measure the work done by the structure on the surrounding fluid.
 
-The accumulated energy can be retrieved using [`calculated_energy`](@ref).
+- **Prescribed/clamped motion energy** (default) -- monitor only the clamped particles by
+    leaving `eachparticle` at its default range
+    `(n_integrated_particles(system) + 1):nparticles(system)`.
+- **Fluid load measurement** -- set `eachparticle=eachparticle(system)` together with
+    `only_compute_force_on_fluid=true` to accumulate the work that the entire structure
+    exerts on the surrounding fluid (useful for drag or lift estimates).
 
-!!! info "TLSPH System Configuration"
-    The [`TotalLagrangianSPHSystem`](@ref) must be created with
-    `use_with_energy_calculator_callback=true` to prepare the system for energy
-    calculation by allocating the necessary cache for clamped particles.
+Internally the callback integrates the instantaneous power, i.e. the dot product between
+the force exerted by the particle and its prescribed velocity, using an explicit Euler
+time integration scheme.
+
+The accumulated value can be retrieved via [`calculated_energy`](@ref).
 
 !!! warning "Experimental implementation"
     This is an experimental feature and may change in future releases.
 
 # Arguments
-- `ELTYPE`: The floating point type used for time and energy calculations.
-            This should match the floating point type used in the
-            [`TotalLagrangianSPHSystem`](@ref), which can be obtained via `eltype(system)`.
+- `system`: The [`TotalLagrangianSPHSystem`](@ref) whose particles should be monitored.
+- `semi`: The [`Semidiscretization`](@ref) that contains `system`.
 
 # Keywords
 - `interval=1`: Interval (in number of time steps) at which to compute the energy.
-                It is recommended to set this either to `1` (every time step) or to a small
-                integer (e.g., `2` or `5`) to reduce time integration errors
-                in the energy calculation.
+                It is recommended to keep this at `1` (every time step) or small (≤ 5)
+                to limit time integration errors in the energy integral.
+- `eachparticle=(n_integrated_particles(system) + 1):nparticles(system)`: Iterator
+                selecting which particles contribute. The default includes all clamped
+                particles in the system; pass `eachparticle(system)` to include every particle.
+- `only_compute_force_on_fluid=false`: When `true`, only interactions with
+                fluid systems are accounted for. Combined with
+                `eachparticle=eachparticle(system)`, this accumulates the work that the
+                entire structure exerts on the fluid, which is useful for drag or lift
+                estimates.
 
 # Examples
-```jldoctest; output = false, setup = :(structure = RectangularShape(0.1, (3, 4), (0.1, 0.0), density=1.0); smoothing_kernel = WendlandC2Kernel{2}(); smoothing_length = 1.0; young_modulus = 1e6; poisson_ratio = 0.3; n_clamped_particles = 0; clamped_particles_motion = nothing)
-# Create TLSPH system with `use_with_energy_calculator_callback=true`
-system = TotalLagrangianSPHSystem(structure, smoothing_kernel, smoothing_length,
-                                  young_modulus, poisson_ratio,
-                                  n_clamped_particles=n_clamped_particles,
-                                  clamped_particles_motion=clamped_particles_motion,
-                                  use_with_energy_calculator_callback=true) # Important!
-
-# Create an energy calculator that runs every 2 time steps
-energy_cb = EnergyCalculatorCallback{eltype(system)}(; interval=2)
+```jldoctest; output = false, setup = :(system = RectangularShape(0.1, (3, 4), (0.1, 0.0), density=1.0); semi = (; systems=(system,), parallelization_backend=SerialBackend()))
+# Create an energy calculator callback that is called every 2 time steps
+energy_cb = EnergyCalculatorCallback(system, semi; interval=2)
 
 # After the simulation, retrieve the calculated energy
 total_energy = calculated_energy(energy_cb)
@@ -60,9 +65,10 @@ struct EnergyCalculatorCallback{T, DV, EP}
     only_compute_force_on_fluid   :: Bool
 end
 
-function EnergyCalculatorCallback{ELTYPE}(system, semi; interval=1,
-                                          eachparticle=eachparticle(system),
-                                          only_compute_force_on_fluid=false) where {ELTYPE <: Real}
+function EnergyCalculatorCallback(system, semi; interval=1,
+                                  eachparticle=(n_integrated_particles(system) + 1):nparticles(system),
+                                  only_compute_force_on_fluid=false)
+    ELTYPE = eltype(system)
     system_index = system_indices(system, semi)
 
     # Allocate buffer to write accelerations for all particles (including clamped ones)
@@ -85,11 +91,11 @@ Get the current calculated energy from an [`EnergyCalculatorCallback`](@ref).
 - `cb`: The `DiscreteCallback` returned by [`EnergyCalculatorCallback`](@ref).
 
 # Examples
-```jldoctest; output = false
-# Before the simulation: create the callback
-energy_cb = EnergyCalculatorCallback{Float64}()
+```jldoctest; output = false, setup = :(system = RectangularShape(0.1, (3, 4), (0.1, 0.0), density=1.0); semi = (; systems=(system,), parallelization_backend=SerialBackend()))
+# Create an energy calculator callback
+energy_cb = EnergyCalculatorCallback(system, semi)
 
-# After the simulation: get the calculated energy
+# After the simulation, retrieve the calculated energy
 total_energy = calculated_energy(energy_cb)
 
 # output