diff --git a/src/callbacks_step/analysis.jl b/src/callbacks_step/analysis.jl
index 414de001..c0eee27a 100644
--- a/src/callbacks_step/analysis.jl
+++ b/src/callbacks_step/analysis.jl
@@ -298,7 +298,7 @@ function analyze_integrals!(io, current_integrals, i, analysis_integrals::NTuple
     quantity = first(analysis_integrals)
     remaining_quantities = Base.tail(analysis_integrals)
 
-    res = analyze(quantity, q, t, semi)
+    res = analyze!(semi, quantity, q, t)
     current_integrals[i] = res
     @printf(io, " %-12s:", pretty_form_utf(quantity))
     @printf(io, "  % 10.8e", res)
@@ -326,17 +326,8 @@ function (cb::DiscreteCallback{Condition, Affect!})(sol) where {Condition,
     return (; l2 = l2_error, linf = linf_error)
 end
 
-function analyze(quantity, q, t, semi::Semidiscretization)
-    integrate_quantity(q -> quantity(q, semi.equations), q, semi)
-end
-
-# The modified entropy/energy of the Svärd-Kalisch and
-# Serre-Green-Naghdi equations
-# takes the whole `q` for every point in space since it requires
-# the derivative of the velocity `v_x`.
-function analyze(quantity::Union{typeof(energy_total_modified), typeof(entropy_modified)},
-                 q, t, semi::Semidiscretization)
-    integrate_quantity(quantity, q, semi)
+function analyze!(semi::Semidiscretization, quantity, q, t)
+    integrate_quantity!(semi.cache.tmp1, quantity, q, semi)
 end
 
 pretty_form_utf(::typeof(waterheight_total)) = "∫η"
diff --git a/src/callbacks_step/relaxation.jl b/src/callbacks_step/relaxation.jl
index 2562e62c..8f85db53 100644
--- a/src/callbacks_step/relaxation.jl
+++ b/src/callbacks_step/relaxation.jl
@@ -69,15 +69,7 @@ end
 
     function relaxation_functional(q, semi)
         @unpack tmp1 = semi.cache
-        # modified entropy from Svärd-Kalisch equations need to take the whole vector `q` for every point in space
-        if relaxation_callback.invariant isa
-           Union{typeof(energy_total_modified), typeof(entropy_modified)}
-            return integrate_quantity!(tmp1, relaxation_callback.invariant, q, semi)
-        else
-            return integrate_quantity!(tmp1,
-                                       q -> relaxation_callback.invariant(q, semi.equations),
-                                       q, semi)
-        end
+        return integrate_quantity!(tmp1, relaxation_callback.invariant, q, semi)
     end
 
     function convex_combination(gamma, old, new)
diff --git a/src/equations/equations.jl b/src/equations/equations.jl
index 92d66e2a..5be8aba9 100644
--- a/src/equations/equations.jl
+++ b/src/equations/equations.jl
@@ -189,7 +189,7 @@ of the correct length `nvariables(equations)`.
 function bathymetry end
 
 """
-        lake_at_rest_error(q, equations::AbstractShallowWaterEquations)
+    lake_at_rest_error(q, equations::AbstractShallowWaterEquations)
 
 Calculate the error for the "lake-at-rest" test case where the
 [`waterheight_total`](@ref) ``\\eta = h + b`` should
@@ -270,12 +270,25 @@ contributions.
 `q_global` is a vector of the primitive variables at ALL nodes.
 `cache` needs to hold the SBP operators used by the `solver` if non-hydrostatic
 terms are present.
+
+Internally, this function allocates a vector for the output and
+calls [`DispersiveShallowWater.energy_total_modified!`](@ref).
 """
 function energy_total_modified(q_global, equations::AbstractShallowWaterEquations, cache)
+    e = similar(q_global.x[begin])
+    return energy_total_modified!(e, q_global, equations, cache)
+end
+
+"""
+    energy_total_modified!(e, q_global, equations::AbstractShallowWaterEquations, cache)
+
+In-place version of [`energy_total_modified`](@ref).
+"""
+function energy_total_modified!(e, q_global, equations::AbstractShallowWaterEquations,
+                                cache)
     # `q_global` is an `ArrayPartition` of the primitive variables at all nodes
     @assert nvariables(equations) == length(q_global.x)
 
-    e = similar(q_global.x[begin])
     for i in eachindex(q_global.x[begin])
         e[i] = energy_total(get_node_vars(q_global, equations, i), equations)
     end
@@ -291,9 +304,20 @@ varnames(::typeof(energy_total_modified), equations) = ("e_modified",)
 Alias for [`energy_total_modified`](@ref).
 """
 @inline function entropy_modified(q_global, equations::AbstractShallowWaterEquations, cache)
-    energy_total_modified(q_global, equations, cache)
+    e = similar(q_global.x[begin])
+    return entropy_modified!(e, q_global, equations, cache)
 end
 
+"""
+    entropy_modified!(e, q_global, equations::AbstractShallowWaterEquations, cache)
+
+In-place version of [`entropy_modified`](@ref).
+"""
+@inline entropy_modified!(e, q_global, equations, cache) = energy_total_modified!(e,
+                                                                                  q_global,
+                                                                                  equations,
+                                                                                  cache)
+
 varnames(::typeof(entropy_modified), equations) = ("U_modified",)
 
 # Add methods to show some information on systems of equations.
@@ -448,3 +472,27 @@ abstract type AbstractSerreGreenNaghdiEquations{NDIMS, NVARS} <:
 const AbstractSerreGreenNaghdiEquations1D = AbstractSerreGreenNaghdiEquations{1}
 include("serre_green_naghdi_1d.jl")
 include("hyperbolic_serre_green_naghdi_1d.jl")
+
+function solve_system_matrix!(dv, system_matrix, rhs,
+                              ::Union{SvaerdKalischEquations1D,
+                                      SerreGreenNaghdiEquations1D},
+                              D1, cache)
+    if issparse(system_matrix)
+        (; factorization) = cache
+        cholesky!(factorization, system_matrix; check = false)
+        if issuccess(factorization)
+            scale_by_mass_matrix!(rhs, D1)
+            # see https://github.com/JoshuaLampert/DispersiveShallowWater.jl/issues/122
+            dv .= factorization \ rhs
+        else
+            # The factorization may fail if the time step is too large
+            # and h becomes negative.
+            fill!(dv, NaN)
+        end
+    else
+        factorization = cholesky!(system_matrix)
+        scale_by_mass_matrix!(rhs, D1)
+        ldiv!(dv, factorization, rhs)
+    end
+    return nothing
+end
diff --git a/src/equations/hyperbolic_serre_green_naghdi_1d.jl b/src/equations/hyperbolic_serre_green_naghdi_1d.jl
index ec8b1d3e..4f90977b 100644
--- a/src/equations/hyperbolic_serre_green_naghdi_1d.jl
+++ b/src/equations/hyperbolic_serre_green_naghdi_1d.jl
@@ -433,9 +433,9 @@ end
 
 # The entropy/energy takes the whole `q` for every point in space
 """
-    energy_total_modified(q_global, equations::HyperbolicSerreGreenNaghdiEquations1D, cache)
+    DispersiveShallowWater.energy_total_modified!(e, q_global, equations::HyperbolicSerreGreenNaghdiEquations1D, cache)
 
-Return the modified total energy of the primitive variables `q_global` for the
+Return the modified total energy `e` of the primitive variables `q_global` for the
 [`HyperbolicSerreGreenNaghdiEquations1D`](@ref).
 It contains additional terms compared to the usual [`energy_total`](@ref)
 modeling non-hydrostatic contributions. The `energy_total_modified`
@@ -449,10 +449,12 @@ the total modified energy is given by
 ```
 
 `q_global` is a vector of the primitive variables at ALL nodes.
+
+See also [`energy_total_modified`](@ref).
 """
-function energy_total_modified(q_global,
-                               equations::HyperbolicSerreGreenNaghdiEquations1D,
-                               cache)
+function energy_total_modified!(e, q_global,
+                                equations::HyperbolicSerreGreenNaghdiEquations1D,
+                                cache)
     # unpack physical parameters and SBP operator `D1`
     g = gravity_constant(equations)
     (; lambda) = equations
@@ -464,8 +466,6 @@ function energy_total_modified(q_global,
     @.. b = equations.eta0 - D
     @.. h = eta - b
 
-    e = zero(h)
-
     # 1/2 g eta^2 + 1/2 h v^2 + 1/6 h^3 w^2 + λ/6 h (1 - H/h)^2
     @.. e = 1 / 2 * g * eta^2 + 1 / 2 * h * v^2 + 1 / 6 * h * w^2 +
             lambda / 6 * h * (1 - H / h)^2
diff --git a/src/equations/serre_green_naghdi_1d.jl b/src/equations/serre_green_naghdi_1d.jl
index a9db959d..c9813f77 100644
--- a/src/equations/serre_green_naghdi_1d.jl
+++ b/src/equations/serre_green_naghdi_1d.jl
@@ -43,7 +43,7 @@ approximation, the Serre-Green-Naghdi equations are
     + \frac{1}{2} g (h^2)_x + g h b_x + \frac{1}{2} h (v^2)_x
     + p_x + \frac{3}{2} \frac{p}{h} b_x + \psi b_x &= 0,\\
   p &= \frac{1}{3} h^3 v_{x}^2 - \frac{1}{3} h^3 v v_{xx}
-    + \frac{1}{2} h^2 v (b_x v)_x,
+    + \frac{1}{2} h^2 v (b_x v)_x,\\
   \psi &= \frac{1}{4} h v (b_x v)_x.
 \end{aligned}
 ```
@@ -57,8 +57,8 @@ References for the Serre-Green-Naghdi system can be found in
   [DOI: 10.1017/S0022112076002425](https://doi.org/10.1017/S0022112076002425)
 
 The semidiscretization implemented here conserves
-- the total water mass (integral of h) as a linear invariant
-- the total momentum (integral of h v) as a nonlinear invariant if the bathymetry is constant
+- the total water mass (integral of ``h``) as a linear invariant
+- the total momentum (integral of ``h v``) as a nonlinear invariant if the bathymetry is constant
 - the total modified energy
 
 for periodic boundary conditions (see Ranocha and Ricchiuto (2024)).
@@ -310,7 +310,7 @@ function create_cache(mesh,
 end
 
 function assemble_system_matrix!(cache, h, D1, D1mat,
-                                 equations::SerreGreenNaghdiEquations1D{BathymetryFlat})
+                                 ::SerreGreenNaghdiEquations1D{BathymetryFlat})
     (; M_h, M_h3_3) = cache
 
     @.. M_h = h
@@ -355,28 +355,6 @@ function assemble_system_matrix!(cache, h, b_x, D1, D1mat,
                      Diagonal(M_h2_bx) * D1mat)
 end
 
-function solve_system_matrix!(dv, system_matrix, rhs,
-                              ::SerreGreenNaghdiEquations1D,
-                              D1, cache)
-    if issparse(system_matrix)
-        (; factorization) = cache
-        cholesky!(factorization, system_matrix; check = false)
-        if issuccess(factorization)
-            scale_by_mass_matrix!(rhs, D1)
-            dv .= factorization \ rhs
-        else
-            # The factorization may fail if the time step is too large
-            # and h becomes negative.
-            fill!(dv, NaN)
-        end
-    else
-        factorization = cholesky!(system_matrix)
-        scale_by_mass_matrix!(rhs, D1)
-        ldiv!(dv, factorization, rhs)
-    end
-    return nothing
-end
-
 # Discretization that conserves
 # - the total water mass (integral of h) as a linear invariant
 # - the total momentum (integral of h v) as a nonlinear invariant for flat bathymetry
@@ -417,7 +395,7 @@ function rhs_sgn_central!(dq, q, equations, source_terms, cache, ::BathymetryFla
     @trixi_timeit timer() "hyperbolic terms" begin
         # Compute all derivatives required below
         (; h, b, h_x, v_x, h2_x, hv_x, v2_x, h2_v_vx_x,
-        h_vx_x, p_x, tmp, M_h, M_h3_3) = cache
+        h_vx_x, p_x, tmp) = cache
 
         @.. b = equations.eta0 - D
         @.. h = eta - b
@@ -499,8 +477,7 @@ function rhs_sgn_upwind!(dq, q, equations, source_terms, cache, ::BathymetryFlat
     @trixi_timeit timer() "hyperbolic terms" begin
         # Compute all derivatives required below
         (; h, b, h_x, v_x, v_x_upwind, h2_x, hv_x, v2_x,
-        h2_v_vx_x, h_vx_x, p_x, tmp,
-        M_h, M_h3_3) = cache
+        h2_v_vx_x, h_vx_x, p_x, tmp) = cache
 
         @.. b = equations.eta0 - D
         @.. h = eta - b
@@ -584,8 +561,7 @@ function rhs_sgn_central!(dq, q, equations, source_terms, cache,
     @trixi_timeit timer() "hyperbolic terms" begin
         # Compute all derivatives required below
         (; h, h_x, v_x, h_hpb_x, b, b_x, hv_x, v2_x,
-        h2_v_vx_x, h_vx_x, p_h, p_x, tmp,
-        M_h_p_h_bx2, M_h3_3, M_h2_bx) = cache
+        h2_v_vx_x, h_vx_x, p_h, p_x, tmp) = cache
         if equations.bathymetry_type isa BathymetryVariable
             (; psi) = cache
         end
@@ -692,8 +668,7 @@ function rhs_sgn_upwind!(dq, q, equations, source_terms, cache,
     @trixi_timeit timer() "hyperbolic terms" begin
         # Compute all derivatives required below
         (; h, h_x, v_x, v_x_upwind, h_hpb_x, b, b_x, hv_x, v2_x,
-        h2_v_vx_x, h_vx_x, p_h, p_0, p_x, tmp,
-        M_h_p_h_bx2, M_h3_3, M_h2_bx) = cache
+        h2_v_vx_x, h_vx_x, p_h, p_0, p_x, tmp) = cache
         if equations.bathymetry_type isa BathymetryVariable
             (; psi) = cache
         end
@@ -806,11 +781,11 @@ end
     return SVector(eta, v, D)
 end
 
-@inline function waterheight_total(q, equations::AbstractSerreGreenNaghdiEquations1D)
+@inline function waterheight_total(q, ::AbstractSerreGreenNaghdiEquations1D)
     return q[1]
 end
 
-@inline function velocity(q, equations::AbstractSerreGreenNaghdiEquations1D)
+@inline function velocity(q, ::AbstractSerreGreenNaghdiEquations1D)
     return q[2]
 end
 
@@ -819,11 +794,11 @@ end
     return equations.eta0 - D
 end
 
-# The entropy/energy takes the whole `q` for every point in space
+# The modified entropy/energy takes the whole `q` for every point in space
 """
-    energy_total_modified(q_global, equations::SerreGreenNaghdiEquations1D, cache)
+    DispersiveShallowWater.energy_total_modified!(e, q_global, equations::SerreGreenNaghdiEquations1D, cache)
 
-Return the modified total energy of the primitive variables `q_global` for the
+Return the modified total energy `e` of the primitive variables `q_global` for the
 [`SerreGreenNaghdiEquations1D`](@ref).
 It contains an additional term containing a
 derivative compared to the usual [`energy_total`](@ref) modeling
@@ -845,24 +820,21 @@ For a [`bathymetry_variable`](@ref) the total modified energy has the additional
 
 `q_global` is a vector of the primitive variables at ALL nodes.
 `cache` needs to hold the SBP operators used by the `solver`.
+
+See also [`energy_total_modified`](@ref).
 """
-function energy_total_modified(q_global,
-                               equations::SerreGreenNaghdiEquations1D,
-                               cache)
+function energy_total_modified!(e, q_global,
+                                equations::SerreGreenNaghdiEquations1D,
+                                cache)
     # unpack physical parameters and SBP operator `D1`
     g = gravity_constant(equations)
     (; D1, h, b, v_x) = cache
 
     # `q_global` is an `ArrayPartition`. It collects the individual arrays for
     # the total water height `eta = h + b` and the velocity `v`.
-    eta, v = q_global.x
-    let D = q_global.x[3]
-        @.. b = equations.eta0 - D
-        @.. h = eta - b
-    end
-
-    N = length(v)
-    e = zeros(eltype(q_global), N)
+    eta, v, D = q_global.x
+    @.. b = equations.eta0 - D
+    @.. h = eta - b
 
     # 1/2 g eta^2 + 1/2 h v^2 + 1/6 h^3 w^2
     # and + 1/8 h (v b_x)^2 for full bathymetry without mild-slope approximation
diff --git a/src/equations/svaerd_kalisch_1d.jl b/src/equations/svaerd_kalisch_1d.jl
index bd8f2793..146e8db1 100644
--- a/src/equations/svaerd_kalisch_1d.jl
+++ b/src/equations/svaerd_kalisch_1d.jl
@@ -1,17 +1,19 @@
 @doc raw"""
-    SvaerdKalischEquations1D(; gravity_constant, eta0 = 0.0,
-                               alpha = 0.0,
-                               beta = 0.2308939393939394,
-                               gamma = 0.04034343434343434)
+    SvaerdKalischEquations1D(bathymetry_type = bathymetry_variable;
+                             gravity_constant, eta0 = 0.0, alpha = 0.0,
+                             beta = 0.2308939393939394, gamma = 0.04034343434343434)
 
-Dispersive system by Svärd and Kalisch in one spatial dimension with spatially varying bathymetry. The equations are given in conservative variables by
+Dispersive system by Svärd and Kalisch in one spatial dimension. The equations for variable bathymetry
+are given in conservative variables by
 ```math
 \begin{aligned}
   h_t + (hv)_x &= (\hat\alpha(\hat\alpha(h + b)_x)_x)_x,\\
   (hv)_t + (hv^2)_x + gh(h + b)_x &= (\hat\alpha v(\hat\alpha(h + b)_x)_x)_x + (\hat\beta v_x)_{xt} + \frac{1}{2}(\hat\gamma v_x)_{xx} + \frac{1}{2}(\hat\gamma v_{xx})_x,
 \end{aligned}
 ```
-where ``\hat\alpha^2 = \alpha\sqrt{gD}D^2``, ``\hat\beta = \beta D^3``, ``\hat\gamma = \gamma\sqrt{gD}D^3``. The coefficients ``\alpha``, ``\beta`` and ``\gamma`` are provided in dimensionless form and ``D = \eta_0 - b`` is the still-water depth and `eta0` is the still-water surface (lake-at-rest).
+where ``\hat\alpha^2 = \alpha\sqrt{gD}D^2``, ``\hat\beta = \beta D^3``, ``\hat\gamma = \gamma\sqrt{gD}D^3``.
+The coefficients ``\alpha``, ``\beta`` and ``\gamma`` are provided in dimensionless form and ``D = \eta_0 - b`` is the
+still-water depth and `eta0` is the still-water surface (lake-at-rest).
 The equations can be rewritten in primitive variables as
 ```math
 \begin{aligned}
@@ -20,14 +22,22 @@ The equations can be rewritten in primitive variables as
 \end{aligned}
 ```
 The unknown quantities of the Svärd-Kalisch equations are the total water height ``\eta`` and the velocity ``v``.
-The gravitational constant is denoted by `g` and the bottom topography (bathymetry) ``b = \eta_0 - D``. The water height above the bathymetry is therefore given by
+The gravitational constant is denoted by `g` and the bottom topography (bathymetry) ``b = \eta_0 - D``.
+The water height above the bathymetry is therefore given by
 ``h = \eta - \eta_0 + D``.
 
+Currently, the equations only support a general variable bathymetry, see [`bathymetry_variable`](@ref).
+
 `SvärdKalischEquations1D` is an alias for `SvaerdKalischEquations1D`.
 
 The equations by Svärd and Kalisch are presented and analyzed in Svärd and Kalisch (2023).
-The semidiscretization implemented here conserves the mass and the energy, is well-balanced for the lake-at-rest state,
-and is developed in Lampert and Ranocha (2024).
+The semidiscretization implemented here conserves
+- the total water (integral of ``h``) as a linear invariant
+- the total momentum (integral of ``h v``) as a nonlinear invariant for flat bathymetry
+- the total modified energy
+
+for periodic boundary conditions (see Lampert, Ranocha).
+Additionally, it is well-balanced for the lake-at-rest stationary solution, see Lampert and Ranocha (2024).
 
 - Magnus Svärd, Henrik Kalisch (2023)
   A novel energy-bounded Boussinesq model and a well-balanced and stable numerical discretization
@@ -36,7 +46,9 @@ and is developed in Lampert and Ranocha (2024).
   Structure-Preserving Numerical Methods for Two Nonlinear Systems of Dispersive Wave Equations
   [DOI: 10.48550/arXiv.2402.16669](https://doi.org/10.48550/arXiv.2402.16669)
 """
-struct SvaerdKalischEquations1D{RealT <: Real} <: AbstractSvaerdKalischEquations{1, 3}
+struct SvaerdKalischEquations1D{Bathymetry <: AbstractBathymetry, RealT <: Real} <:
+       AbstractSvaerdKalischEquations{1, 3}
+    bathymetry_type::Bathymetry # type of bathymetry
     gravity::RealT # gravitational constant
     eta0::RealT    # constant still-water surface
     alpha::RealT   # coefficient
@@ -46,9 +58,10 @@ end
 
 const SvärdKalischEquations1D = SvaerdKalischEquations1D
 
-function SvaerdKalischEquations1D(; gravity_constant, eta0 = 0.0, alpha = 0.0,
+function SvaerdKalischEquations1D(bathymetry_type = bathymetry_variable;
+                                  gravity_constant, eta0 = 0.0, alpha = 0.0,
                                   beta = 0.2308939393939394, gamma = 0.04034343434343434)
-    SvaerdKalischEquations1D(gravity_constant, eta0, alpha, beta, gamma)
+    SvaerdKalischEquations1D(bathymetry_type, gravity_constant, eta0, alpha, beta, gamma)
 end
 
 varnames(::typeof(prim2prim), ::SvaerdKalischEquations1D) = ("η", "v", "D")
@@ -70,16 +83,17 @@ References:
   [link](https://repository.tudelft.nl/islandora/object/uuid:c2091d53-f455-48af-a84b-ac86680455e9/datastream/OBJ/download)
 """
 function initial_condition_dingemans(x, t, equations::SvaerdKalischEquations1D, mesh)
+    g = gravity_constant(equations)
     h0 = 0.8
     A = 0.02
     # omega = 2*pi/(2.02*sqrt(2))
-    k = 0.8406220896381442 # precomputed result of find_zero(k -> omega^2 - equations.gravity * k * tanh(k * h0), 1.0) using Roots.jl
+    k = 0.8406220896381442 # precomputed result of find_zero(k -> omega^2 - g * k * tanh(k * h0), 1.0) using Roots.jl
     if x < -30.5 * pi / k || x > -8.5 * pi / k
         h = 0.0
     else
         h = A * cos(k * x)
     end
-    v = sqrt(equations.gravity / k * tanh(k * h0)) * h / h0
+    v = sqrt(g / k * tanh(k * h0)) * h / h0
     if 11.01 <= x && x < 23.04
         b = 0.6 * (x - 11.01) / (23.04 - 11.01)
     elseif 23.04 <= x && x < 27.04
@@ -115,8 +129,8 @@ end
 A smooth manufactured solution in combination with [`initial_condition_manufactured`](@ref).
 """
 function source_terms_manufactured(q, x, t, equations::SvaerdKalischEquations1D)
-    g = equations.gravity
-    eta0 = equations.eta0
+    g = gravity_constant(equations)
+    eta0 = still_water_surface(q, equations)
     alpha = equations.alpha
     beta = equations.beta
     gamma = equations.gamma
@@ -174,129 +188,188 @@ function create_cache(mesh, equations::SvaerdKalischEquations1D,
                       solver, initial_condition,
                       ::BoundaryConditionPeriodic,
                       RealT, uEltype)
+    D1 = solver.D1
     #  Assume D is independent of time and compute D evaluated at mesh points once.
     D = Array{RealT}(undef, nnodes(mesh))
     x = grid(solver)
     for i in eachnode(solver)
         D[i] = still_waterdepth(initial_condition(x[i], 0.0, equations, mesh), equations)
     end
-    h = Array{RealT}(undef, nnodes(mesh))
-    hv = similar(h)
-    alpha_hat = sqrt.(equations.alpha * sqrt.(equations.gravity * D) .* D .^ 2)
+    g = gravity_constant(equations)
+    h = ones(RealT, nnodes(mesh))
+    hv = zero(h)
+    b = zero(h)
+    eta_x = zero(h)
+    v_x = zero(h)
+    alpha_eta_x_x = zero(h)
+    y_x = zero(h)
+    v_y_x = zero(h)
+    yv_x = zero(h)
+    vy = zero(h)
+    vy_x = zero(h)
+    y_v_x = zero(h)
+    h_v_x = zero(h)
+    hv2_x = zero(h)
+    v_xx = zero(h)
+    gamma_v_xx_x = zero(h)
+    gamma_v_x_xx = zero(h)
+    alpha_hat = sqrt.(equations.alpha * sqrt.(g * D) .* D .^ 2)
     beta_hat = equations.beta * D .^ 3
-    gamma_hat = equations.gamma * sqrt.(equations.gravity * D) .* D .^ 3
-    tmp2 = similar(h)
-    M = mass_matrix(solver.D1)
-    if solver.D1 isa PeriodicDerivativeOperator ||
-       solver.D1 isa UniformPeriodicCoupledOperator
-        D1_central = solver.D1
-        sparse_D1 = sparse(D1_central)
-        # We use the periodic SBP property
-        #   M * sparse_D1 == -sparse_D1' * M
-        # to avoid possible floating point errors in the symmetry of the matrix.
-        minus_MD1betaD1 = sparse_D1' * M * Diagonal(beta_hat) * sparse_D1
-    elseif solver.D1 isa PeriodicUpwindOperators
-        D1_central = solver.D1.central
-        # We use the periodic upwind SBP property
-        #   M * sparse_D1plus == -sparse_D1minus' * M
-        # to avoid possible floating point errors in the symmetry of the matrix.
-        sparse_D1minus = sparse(solver.D1.minus)
-        minus_MD1betaD1 = sparse_D1minus' * M * Diagonal(beta_hat) * sparse_D1minus
+    gamma_hat = equations.gamma * sqrt.(g * D) .* D .^ 3
+    tmp2 = zero(h)
+    M = mass_matrix(D1)
+    M_h = zero(h)
+    M_beta = copy(beta_hat)
+    scale_by_mass_matrix!(M_beta, D1)
+    if D1 isa PeriodicDerivativeOperator ||
+       D1 isa UniformPeriodicCoupledOperator
+        D1_central = D1
+        D1mat = sparse(D1_central)
+        minus_MD1betaD1 = D1mat' * Diagonal(M_beta) * D1mat
+        system_matrix = let cache = (; M_h, minus_MD1betaD1)
+            assemble_system_matrix!(cache, h,
+                                    D1, D1mat, equations)
+        end
+    elseif D1 isa PeriodicUpwindOperators
+        D1_central = D1.central
+        D1mat_minus = sparse(D1.minus)
+        minus_MD1betaD1 = D1mat_minus' * Diagonal(M_beta) * D1mat_minus
+        system_matrix = let cache = (; M_h, minus_MD1betaD1)
+            assemble_system_matrix!(cache, h,
+                                    D1, D1mat_minus, equations)
+        end
     else
-        @error "unknown type of first-derivative operator: $(typeof(solver.D1))"
+        @error "unknown type of first-derivative operator: $(typeof(D1))"
+    end
+    factorization = cholesky(system_matrix)
+    cache = (; factorization, minus_MD1betaD1, D, h, hv, b, eta_x, v_x,
+             alpha_eta_x_x, y_x, v_y_x, yv_x, vy, vy_x, y_v_x, h_v_x, hv2_x, v_xx,
+             gamma_v_xx_x, gamma_v_x_xx,
+             alpha_hat, beta_hat, gamma_hat, tmp2, D1_central, M, D1, M_h)
+    if D1 isa PeriodicUpwindOperators
+        eta_x_upwind = zero(h)
+        v_x_upwind = zero(h)
+        cache = (; cache..., eta_x_upwind, v_x_upwind)
     end
-    factorization = cholesky(Symmetric(M * Diagonal(ones(nnodes(mesh))) + minus_MD1betaD1))
-    return (factorization = factorization, minus_MD1betaD1 = minus_MD1betaD1, D = D, h = h,
-            hv = hv, alpha_hat = alpha_hat, gamma_hat = gamma_hat,
-            tmp2 = tmp2, D1_central = D1_central, M = M, D1 = solver.D1)
+    return cache
+end
+
+function assemble_system_matrix!(cache, h, D1, D1mat,
+                                 ::SvaerdKalischEquations1D)
+    (; M_h, minus_MD1betaD1) = cache
+
+    @.. M_h = h
+    scale_by_mass_matrix!(M_h, D1)
+
+    # Floating point errors accumulate a bit and the system matrix
+    # is not necessarily perfectly symmetric but only up to
+    # round-off errors. We wrap it here to avoid issues with the
+    # factorization.
+    return Symmetric(Diagonal(M_h) + minus_MD1betaD1)
 end
 
-# Discretization that conserves the mass (for eta and for flat bottom hv) and the energy for periodic boundary conditions, see
+# Discretization that conserves
+# - the total water (integral of h) as a linear invariant
+# - the total momentum (integral of h v) as a nonlinear invariant for flat bathymetry
+# - the total modified energy
+# for periodic boundary conditions, see
 # - Joshua Lampert and Hendrik Ranocha (2024)
 #   Structure-Preserving Numerical Methods for Two Nonlinear Systems of Dispersive Wave Equations
 #   [DOI: 10.48550/arXiv.2402.16669](https://doi.org/10.48550/arXiv.2402.16669)
+# TODO: Simplify for the case of flat bathymetry and use higher-order operators
 function rhs!(dq, q, t, mesh, equations::SvaerdKalischEquations1D,
               initial_condition, ::BoundaryConditionPeriodic, source_terms,
               solver, cache)
-    @unpack factorization, minus_MD1betaD1, D, h, hv, alpha_hat, gamma_hat, tmp1, tmp2, D1_central, M = cache
+    (; D, h, hv, b, eta_x, v_x, alpha_eta_x_x, y_x, v_y_x, yv_x, vy, vy_x,
+    y_v_x, h_v_x, hv2_x, v_xx, gamma_v_xx_x, gamma_v_x_xx, alpha_hat, gamma_hat,
+    tmp1, tmp2, D1_central, M, D1) = cache
 
-    eta = q.x[1]
-    v = q.x[2]
-    deta = dq.x[1]
-    dv = dq.x[2]
-    dD = dq.x[3]
+    g = gravity_constant(equations)
+    eta, v = q.x
+    deta, dv, dD = dq.x
     fill!(dD, zero(eltype(dD)))
 
     @trixi_timeit timer() "deta hyperbolic" begin
-        @. h = eta + D - equations.eta0
-        @. hv = h * v
-
-        if solver.D1 isa PeriodicDerivativeOperator ||
-           solver.D1 isa UniformPeriodicCoupledOperator
-            D1eta = D1_central * eta
-            D1v = D1_central * v
-            tmp1 = alpha_hat .* (D1_central * (alpha_hat .* D1eta))
-            vD1y = v .* (D1_central * tmp1)
-            D1vy = D1_central * (v .* tmp1)
-            yD1v = tmp1 .* D1v
-            @. tmp2 = tmp1 - hv
+        @.. b = equations.eta0 - D
+        @.. h = eta - b
+        @.. hv = h * v
+
+        if D1 isa PeriodicDerivativeOperator ||
+           D1 isa UniformPeriodicCoupledOperator
+            mul!(eta_x, D1_central, eta)
+            mul!(v_x, D1_central, v)
+            @.. tmp1 = alpha_hat * eta_x
+            mul!(alpha_eta_x_x, D1_central, tmp1)
+            @.. tmp1 = alpha_hat * alpha_eta_x_x
+            mul!(y_x, D1_central, tmp1)
+            @.. v_y_x = v * y_x
+            @.. vy = v * tmp1
+            mul!(vy_x, D1_central, vy)
+            @.. y_v_x = tmp1 * v_x
+            @.. tmp2 = tmp1 - hv
             mul!(deta, D1_central, tmp2)
-        elseif solver.D1 isa PeriodicUpwindOperators
-            D1eta = D1_central * eta
-            D1v = D1_central * v
-            tmp1 = alpha_hat .* (solver.D1.minus * (alpha_hat .* (solver.D1.plus * eta)))
-            vD1y = v .* (solver.D1.minus * tmp1)
-            D1vy = solver.D1.minus * (v .* tmp1)
-            yD1v = tmp1 .* (solver.D1.plus * v)
-            deta[:] = solver.D1.minus * tmp1 - D1_central * hv
+        elseif D1 isa PeriodicUpwindOperators
+            @unpack eta_x_upwind, v_x_upwind = cache
+            mul!(eta_x, D1_central, eta)
+            mul!(v_x, D1_central, v)
+            mul!(eta_x_upwind, D1.plus, eta)
+            @.. tmp1 = alpha_hat * eta_x_upwind
+            mul!(alpha_eta_x_x, D1_central, tmp1)
+            @.. tmp1 = alpha_hat * alpha_eta_x_x
+            mul!(y_x, D1.minus, tmp1)
+            @.. v_y_x = v * y_x
+            @.. vy = v * tmp1
+            mul!(vy_x, D1.minus, vy)
+            mul!(v_x_upwind, D1.plus, v)
+            @.. y_v_x = tmp1 * v_x_upwind
+            # deta[:] = D1.minus * tmp1 - D1_central * hv
+            mul!(deta, D1.minus, tmp1)
+            mul!(deta, D1_central, hv, -1.0, 1.0)
         else
-            @error "unknown type of first derivative operator: $(typeof(solver.D1))"
+            @error "unknown type of first derivative operator: $(typeof(D1))"
         end
     end
 
     # split form
     @trixi_timeit timer() "dv hyperbolic" begin
-        D1_hv = D1_central * hv
-        D1_hv2 = D1_central * (hv .* v)
-        D1_gamma_hat_D2_v = D1_central * (gamma_hat .* (solver.D2 * v))
-        D2_gamma_hat_D1_v = solver.D2 * (gamma_hat .* D1v)
-        @. dv = -(0.5 * (D1_hv2 + hv * D1v - v * D1_hv) +
-                  equations.gravity * h * D1eta +
-                  0.5 * (vD1y - D1vy - yD1v) -
-                  0.5 * D1_gamma_hat_D2_v -
-                  0.5 * D2_gamma_hat_D1_v)
+        mul!(h_v_x, D1_central, hv)
+        @.. tmp1 = hv * v
+        mul!(hv2_x, D1_central, tmp1)
+        mul!(v_xx, solver.D2, v)
+        @.. tmp1 = gamma_hat * v_xx
+        mul!(gamma_v_xx_x, D1_central, tmp1)
+        @.. tmp1 = gamma_hat * v_x
+        mul!(gamma_v_x_xx, solver.D2, tmp1)
+        @.. dv = -(0.5 * (hv2_x + hv * v_x - v * h_v_x) +
+                   g * h * eta_x +
+                   0.5 * (v_y_x - vy_x - y_v_x) -
+                   0.5 * gamma_v_xx_x -
+                   0.5 * gamma_v_x_xx)
     end
 
     # no split form
     #     dv[:] = -(D1_central * (hv .* v) - v .* (D1_central * hv)+
-    #               equations.gravity * h .* D1eta +
-    #               vD1y - D1vy -
-    #               0.5 * D1_central * (gamma_hat .* (solver.D2 * v)) -
-    #               0.5 * solver.D2 * (gamma_hat .* D1v))
+    #               g * h .* eta_x +
+    #               vy_x - v_y_c -
+    #               0.5 * gamma_v_xx_x -
+    #               0.5 * gamma_v_x_xx)
 
     @trixi_timeit timer() "source terms" calc_sources!(dq, q, t, source_terms, equations,
                                                        solver)
-    @trixi_timeit timer() "dv elliptic" begin
-        # decompose M * (h - D1betaD1) because it is guaranteed to be symmetric and pos. def.,
-        # while (h - D1betaD1) is not necessarily
-        hmD1betaD1 = Symmetric(M * Diagonal(h) + minus_MD1betaD1)
-        # If the time integration method takes a too large step, h become become negative
-        # letting the factorization fail. In this case the solution is set to `NaN` to force
-        # rejecting the step
-        cholesky!(factorization, hmD1betaD1, check = false)
-        if issuccess(factorization)
-            mul!(tmp1, M, dv)
-            dv[:] = factorization \ tmp1
-        else
-            fill!(dv, NaN)
-        end
+    @trixi_timeit timer() "assemble system matrix" begin
+        system_matrix = assemble_system_matrix!(cache, h, D1, D1_central, equations)
+    end
+    @trixi_timeit timer() "solving elliptic system" begin
+        tmp1 .= dv
+        solve_system_matrix!(dv, system_matrix, tmp1,
+                             equations, D1, cache)
     end
 
     return nothing
 end
 
 @inline function prim2cons(q, equations::SvaerdKalischEquations1D)
-    eta, v, D = q
+    eta, v = q
 
     b = bathymetry(q, equations)
     h = eta - b
@@ -313,11 +386,11 @@ end
     return SVector(eta, v, D)
 end
 
-@inline function waterheight_total(q, equations::SvaerdKalischEquations1D)
+@inline function waterheight_total(q, ::SvaerdKalischEquations1D)
     return q[1]
 end
 
-@inline function velocity(q, equations::SvaerdKalischEquations1D)
+@inline function velocity(q, ::SvaerdKalischEquations1D)
     return q[2]
 end
 
@@ -326,41 +399,44 @@ end
     return equations.eta0 - D
 end
 
-@inline entropy(u, equations::SvaerdKalischEquations1D) = energy_total(u, equations)
-
-# The modified entropy/total energy takes the whole `q` for every point in space
+# The modified entropy/energy takes the whole `q` for every point in space
 """
-    energy_total_modified(q_global, equations::SvaerdKalischEquations1D, cache)
+    DispersiveShallowWater.energy_total_modified!(e, q_global, equations::SvaerdKalischEquations1D, cache)
 
-Return the modified total energy of the primitive variables `q_global` for the
+Return the modified total energy `e` of the primitive variables `q_global` for the
 [`SvaerdKalischEquations1D`](@ref). It contains an additional term containing a
-derivative compared to the usual [`energy_total`](@ref). The `energy_total_modified`
-is a conserved quantity of the Svärd-Kalisch equations given by
+derivative compared to the usual [`energy_total`](@ref) modeling
+non-hydrostatic contributions. The `energy_total_modified`
+is a conserved quantity (for periodic boundary conditions).
+
+It is given by
 ```math
-\\frac{1}{2} g h^2 + \\frac{1}{2} h v^2 + \\frac{1}{2} \\hat\\beta v_x^2.
+\\frac{1}{2} g \\eta^2 + \\frac{1}{2} h v^2 + \\frac{1}{2} \\hat\\beta v_x^2.
 ```
 
 `q_global` is a vector of the primitive variables at ALL nodes.
-`cache` needs to hold the first-derivative SBP operator `D1`.
+`cache` needs to hold the SBP operators used by the `solver`.
+
+See also [`energy_total_modified`](@ref).
 """
-@inline function energy_total_modified(q_global, equations::SvaerdKalischEquations1D, cache)
-    # Need to compute new beta_hat, do not use the old one from the `cache`
-    v = q_global.x[2]
-    D = q_global.x[3]
-    N = length(v)
-    e_modified = zeros(eltype(q_global), N)
-    beta_hat = equations.beta * D .^ 3
-    if cache.D1 isa PeriodicDerivativeOperator ||
-       cache.D1 isa UniformPeriodicCoupledOperator
-        tmp = 0.5 * beta_hat .* ((cache.D1 * v) .^ 2)
-    elseif cache.D1 isa PeriodicUpwindOperators
-        tmp = 0.5 * beta_hat .* ((cache.D1.minus * v) .^ 2)
+@inline function energy_total_modified!(e, q_global, equations::SvaerdKalischEquations1D,
+                                        cache)
+    # unpack physical parameters and SBP operator `D1`
+    g = gravity_constant(equations)
+    (; D1, h, b, v_x, beta_hat) = cache
+
+    # `q_global` is an `ArrayPartition`. It collects the individual arrays for
+    # the total water height `eta = h + b` and the velocity `v`.
+    eta, v, D = q_global.x
+    @.. b = equations.eta0 - D
+    @.. h = eta - b
+
+    if D1 isa PeriodicUpwindOperators
+        mul!(v_x, D1.minus, v)
     else
-        @error "unknown type of first-derivative operator: $(typeof(cache.D1))"
+        mul!(v_x, D1, v)
     end
-    for i in 1:N
-        e_modified[i] = energy_total(get_node_vars(q_global, equations, i), equations) +
-                        tmp[i]
-    end
-    return e_modified
+
+    @.. e = 1 / 2 * g * eta^2 + 1 / 2 * h * v^2 + 1 / 2 * beta_hat * v_x^2
+    return e
 end
diff --git a/src/semidiscretization.jl b/src/semidiscretization.jl
index 60e6fa33..f89a0bd8 100644
--- a/src/semidiscretization.jl
+++ b/src/semidiscretization.jl
@@ -139,26 +139,23 @@ end
 
 function integrate_quantity!(quantity, func, q, semi::Semidiscretization)
     for i in eachnode(semi)
-        quantity[i] = func(get_node_vars(q, semi.equations, i))
+        quantity[i] = func(get_node_vars(q, semi.equations, i), semi.equations)
     end
     integrate(quantity, semi)
 end
 
+# Obtain the function, which has an additional `!` appended to the name
+inplace_version(f) = getfield(@__MODULE__, Symbol(string(nameof(f)) * "!"))
+
 # The entropy/energy of the Svärd-Kalisch and Serre-Green-Naghdi equations
 # takes the whole `q` for every point in space since it requires
 # the derivative of the velocity `v_x`.
-function integrate_quantity(func::Union{typeof(energy_total_modified),
-                                        typeof(entropy_modified)}, q,
-                            semi::Semidiscretization)
-    quantity = func(q, semi.equations, semi.cache)
-    integrate(quantity, semi)
-end
-
 function integrate_quantity!(quantity,
                              func::Union{typeof(energy_total_modified),
                                          typeof(entropy_modified)}, q,
                              semi::Semidiscretization)
-    integrate_quantity(func, q, semi)
+    inplace_version(func)(quantity, q, semi.equations, semi.cache)
+    integrate(quantity, semi)
 end
 
 @inline function mesh_equations_solver(semi::Semidiscretization)
diff --git a/test/test_unit.jl b/test/test_unit.jl
index fa309a65..ca8a38af 100644
--- a/test/test_unit.jl
+++ b/test/test_unit.jl
@@ -114,6 +114,8 @@ using SparseArrays: sparse, SparseMatrixCSC
             prim2cons,
             prim2prim,
             prim2phys,
+            energy_total_modified,
+            entropy_modified,
         ]
         for conversion in conversion_functions
             @test DispersiveShallowWater.varnames(conversion, equations) isa Tuple
@@ -143,10 +145,12 @@ using SparseArrays: sparse, SparseMatrixCSC
             _, _, _, cache = @inferred DispersiveShallowWater.mesh_equations_solver_cache(semi)
             e_modified = @inferred energy_total_modified(q, equations, cache)
             e_modified_total = @inferred DispersiveShallowWater.integrate(e_modified, semi)
-            e_total = @inferred DispersiveShallowWater.integrate_quantity(q -> energy_total(q,
-                                                                                            equations),
+            e_total = @inferred DispersiveShallowWater.integrate_quantity(energy_total,
                                                                           q, semi)
             @test isapprox(e_modified_total, e_total)
+            U_modified = @inferred entropy_modified(q, equations, cache)
+            U_modified_total = @inferred DispersiveShallowWater.integrate(U_modified, semi)
+            @test isapprox(U_modified_total, e_modified_total)
         end
     end
 
@@ -165,6 +169,8 @@ using SparseArrays: sparse, SparseMatrixCSC
             prim2cons,
             prim2prim,
             prim2phys,
+            energy_total_modified,
+            entropy_modified,
         ]
         for conversion in conversion_functions
             @test DispersiveShallowWater.varnames(conversion, equations) isa Tuple
@@ -194,10 +200,12 @@ using SparseArrays: sparse, SparseMatrixCSC
             _, _, _, cache = @inferred DispersiveShallowWater.mesh_equations_solver_cache(semi)
             e_modified = @inferred energy_total_modified(q, equations, cache)
             e_modified_total = @inferred DispersiveShallowWater.integrate(e_modified, semi)
-            e_total = @inferred DispersiveShallowWater.integrate_quantity(q -> energy_total(q,
-                                                                                            equations),
+            e_total = @inferred DispersiveShallowWater.integrate_quantity(energy_total,
                                                                           q, semi)
             @test isapprox(e_modified_total, e_total)
+            U_modified = @inferred entropy_modified(q, equations, cache)
+            U_modified_total = @inferred DispersiveShallowWater.integrate(U_modified, semi)
+            @test isapprox(U_modified_total, e_modified_total)
         end
     end
 
@@ -236,6 +244,27 @@ using SparseArrays: sparse, SparseMatrixCSC
         @test @inferred(still_water_surface(q, equations)) == 0.0
         @test isapprox(@inferred(energy_total(q, equations)), 8740.42)
         @test @inferred(prim2phys(q, equations)) == @inferred(prim2prim(q, equations))
+
+        @testset "energy_total_modified" begin
+            initial_condition = initial_condition_manufactured
+            boundary_conditions = boundary_condition_periodic
+            mesh = @inferred Mesh1D(-1.0, 1.0, 10)
+            solver = Solver(mesh, 4)
+            semi = @inferred Semidiscretization(mesh, equations, initial_condition,
+                                                solver; boundary_conditions)
+            q = @inferred DispersiveShallowWater.compute_coefficients(initial_condition,
+                                                                      0.0, semi)
+            _, _, _, cache = @inferred DispersiveShallowWater.mesh_equations_solver_cache(semi)
+            e_modified = @inferred energy_total_modified(q, equations, cache)
+            e_modified_total = @inferred DispersiveShallowWater.integrate(e_modified, semi)
+            e_total = @inferred DispersiveShallowWater.integrate_quantity(energy_total,
+                                                                          q, semi)
+            @test isapprox(e_modified_total, 1450.0018635214328)
+            @test isapprox(e_total, 7.405000000000001)
+            U_modified = @inferred entropy_modified(q, equations, cache)
+            U_modified_total = @inferred DispersiveShallowWater.integrate(U_modified, semi)
+            @test isapprox(U_modified_total, e_modified_total)
+        end
     end
 
     @testset "SerreGreenNaghdiEquations1D" begin
@@ -253,6 +282,8 @@ using SparseArrays: sparse, SparseMatrixCSC
             prim2cons,
             prim2prim,
             prim2phys,
+            energy_total_modified,
+            entropy_modified,
         ]
         for conversion in conversion_functions
             @test DispersiveShallowWater.varnames(conversion, equations) isa Tuple
@@ -267,6 +298,27 @@ using SparseArrays: sparse, SparseMatrixCSC
         @test @inferred(discharge(q, equations)) == 84.0
         @test @inferred(still_water_surface(q, equations)) == 0.0
         @test @inferred(prim2phys(q, equations)) == @inferred(prim2prim(q, equations))
+
+        @testset "energy_total_modified" begin
+            initial_condition = initial_condition_convergence_test
+            boundary_conditions = boundary_condition_periodic
+            mesh = @inferred Mesh1D(-1.0, 1.0, 10)
+            solver = Solver(mesh, 4)
+            semi = @inferred Semidiscretization(mesh, equations, initial_condition,
+                                                solver; boundary_conditions)
+            q = @inferred DispersiveShallowWater.compute_coefficients(initial_condition,
+                                                                      0.0, semi)
+            _, _, _, cache = @inferred DispersiveShallowWater.mesh_equations_solver_cache(semi)
+            e_modified = @inferred energy_total_modified(q, equations, cache)
+            e_modified_total = @inferred DispersiveShallowWater.integrate(e_modified, semi)
+            e_total = @inferred DispersiveShallowWater.integrate_quantity(energy_total,
+                                                                          q, semi)
+            @test isapprox(e_modified_total, 14.303587674490101)
+            @test isapprox(e_total, 14.301514636021535)
+            U_modified = @inferred entropy_modified(q, equations, cache)
+            U_modified_total = @inferred DispersiveShallowWater.integrate(U_modified, semi)
+            @test isapprox(U_modified_total, e_modified_total)
+        end
     end
 
     @testset "HyperbolicSerreGreenNaghdiEquations1D" begin
@@ -285,6 +337,8 @@ using SparseArrays: sparse, SparseMatrixCSC
             prim2cons,
             prim2prim,
             prim2phys,
+            energy_total_modified,
+            entropy_modified,
         ]
         for conversion in conversion_functions
             @test DispersiveShallowWater.varnames(conversion, equations) isa Tuple
@@ -299,6 +353,27 @@ using SparseArrays: sparse, SparseMatrixCSC
         @test @inferred(discharge(q, equations)) == 84.0
         @test @inferred(still_water_surface(q, equations)) == 0.0
         @test @inferred(prim2phys(q, equations)) == [42.0, 2.0, 0.0]
+
+        @testset "energy_total_modified" begin
+            initial_condition = initial_condition_soliton
+            boundary_conditions = boundary_condition_periodic
+            mesh = @inferred Mesh1D(-1.0, 1.0, 10)
+            solver = Solver(mesh, 4)
+            semi = @inferred Semidiscretization(mesh, equations, initial_condition,
+                                                solver; boundary_conditions)
+            q = @inferred DispersiveShallowWater.compute_coefficients(initial_condition,
+                                                                      0.0, semi)
+            _, _, _, cache = @inferred DispersiveShallowWater.mesh_equations_solver_cache(semi)
+            e_modified = @inferred energy_total_modified(q, equations, cache)
+            e_modified_total = @inferred DispersiveShallowWater.integrate(e_modified, semi)
+            e_total = @inferred DispersiveShallowWater.integrate_quantity(energy_total,
+                                                                          q, semi)
+            @test isapprox(e_modified_total, 14.303814990428117)
+            @test isapprox(e_total, 14.301514636021535)
+            U_modified = @inferred entropy_modified(q, equations, cache)
+            U_modified_total = @inferred DispersiveShallowWater.integrate(U_modified, semi)
+            @test isapprox(U_modified_total, e_modified_total)
+        end
     end
 
     @testset "AnalysisCallback" begin