diff --git a/src/Trixi.jl b/src/Trixi.jl
index 4b7bc79d6c6..2d24b43ddce 100644
--- a/src/Trixi.jl
+++ b/src/Trixi.jl
@@ -244,9 +244,9 @@ export cons2cons, cons2prim, prim2cons, cons2macroscopic, cons2state, cons2mean,
        cons2entropy, entropy2cons
 export density, pressure, density_pressure, velocity, global_mean_vars,
        equilibrium_distribution, waterheight, waterheight_pressure
-export entropy, energy_total, energy_kinetic, energy_internal,
-       energy_magnetic, cross_helicity, magnetic_field, divergence_cleaning_field,
-       enstrophy, vorticity
+export entropy, energy_total, energy_kinetic, energy_internal, energy_magnetic,
+       cross_helicity,
+       enstrophy, magnetic_field, divergence_cleaning_field, enstrophy_multi_euler, vorticity
 export lake_at_rest_error
 export ncomponents, eachcomponent
 
diff --git a/src/callbacks_step/analysis_dg2d.jl b/src/callbacks_step/analysis_dg2d.jl
index fa18c5af63a..6ae9cf65f5b 100644
--- a/src/callbacks_step/analysis_dg2d.jl
+++ b/src/callbacks_step/analysis_dg2d.jl
@@ -246,10 +246,135 @@ function integrate(func::Func, u,
                                UnstructuredMesh2D, P4estMesh{2}, P4estMeshView{2},
                                T8codeMesh{2}},
                    equations, dg::DG, cache; normalize = true) where {Func}
-    integrate_via_indices(u, mesh, equations, dg, cache;
-                          normalize = normalize) do u, i, j, element, equations, dg
-        u_local = get_node_vars(u, equations, dg, i, j, element)
-        return func(u_local, equations)
+    @unpack boundaries = cache
+    m = methods(func)
+    if (m[1].nargs == 2) || (func == cons2cons)
+        return integrate_via_indices(u, mesh, equations, dg, cache;
+                                     normalize = normalize) do u, i, j, element,
+                                                               equations, dg
+            u_local = get_node_vars(u, equations, dg, i, j, element)
+
+            func(u_local, equations)
+        end
+    end
+    if (m[1].nargs == 4) && (func != cons2cons)
+        gradients_x = DGSpaceDerivative_WeakForm!(dg, cache, u, 1, equations, mesh)
+        gradients_y = DGSpaceDerivative_WeakForm!(dg, cache, u, 2, equations, mesh)
+        return integrate_via_indices(u, mesh, equations, dg, cache;
+                                     normalize = normalize) do u, i, j, element,
+                                                               equations, dg
+            u_local = get_node_vars(u, equations, dg, i, j, element)
+            gradients_local = Vector([gradients_x[:, i, j, element], gradients_y[:, i, j, element]])
+
+            func(u_local, gradients_local, equations)
+        end
+    end
+end
+
+# Return the derivatives of the interpolatied polynomial.
+function lagrange_derivatives(x, y)
+    n = length(x)
+    dydx = zeros(n)
+
+    for i in 1:n
+        for j in 1:n
+            if i != j
+                term = y[j] / (x[i] - x[j])
+                for k in 1:n
+                    if k != i && k != j
+                        term *= (x[i] - x[k]) / (x[j] - x[k])
+                    end
+                end
+                dydx[i] += term
+            end
+        end
+    end
+
+    return dydx
+end
+
+# Andrew's functions for computing the derivatives.
+function DGSpaceDerivative_WeakForm!(dg,
+                                     cache,
+                                     u,
+                                     direction::Int,
+                                     equations,
+                                     mesh)
+    # Get the required variables.
+    @unpack derivative_dhat, weights = dg.basis
+
+    # Compute the surface flux terms.
+    surface_flux_values = similar(u)
+    Trixi.calc_surface_integral!(surface_flux_values, u, mesh, equations, dg.surface_integral, dg, cache)
+
+    # Translations.
+    D = derivative_dhat
+    spA_Dhat = D
+
+    n_vars, Np, _, n_elements = size(u)
+    gradients = similar(u)
+
+    u_local = zeros((n_vars, Np, Np))
+
+    if direction == 1
+        for elem in 1:n_elements
+            # We work in primitive variables here.
+            for i in 1:Np
+                for j in 1:Np
+                    u_local[:, i, j] = cons2prim(u[:, i, j, elem], equations)
+                end
+            end
+            # Volume contributions (weak form)
+            for v in 1:n_vars
+                # Interpolate polynomial.
+                for j in 1:Np
+                    x = cache.elements.node_coordinates[1, :, j, elem]
+                    y = u_local[v, :, j]
+                    gradients[v, :, j, elem] .= lagrange_derivatives(x, y)
+                end
+
+#                 # Contract D in the ξ̂ (first spatial) direction
+#                 gradients[v, :, :, elem] .= D * u[v, :, :, elem]
+#                 # Average over element.
+
+                average = sum(gradients[v, :, :, elem]/Np^2)
+                gradients[v, :, :, elem] .= average
+            end
+        end
+    elseif direction == 2
+        for elem in 1:n_elements
+            # We work in primitive variables here.
+            for i in 1:Np
+                for j in 1:Np
+                    u_local[:, i, j] = cons2prim(u[:, i, j, elem], equations)
+                end
+            end
+            # Volume contributions (weak form)
+            for v in 1:n_vars
+                # Interpolate polynomial.
+                for i in 1:Np
+                    x = cache.elements.node_coordinates[2, i, :, elem]
+                    y = u_local[v, i, :]
+                    gradients[v, i, :, elem] .= lagrange_derivatives(x, y)
+                end
+
+#                 # Contract D in the ξ̂ (first spatial) direction
+#                 gradients[v, :, :, elem] .= u[v, :, :, elem] * D'
+
+                # Average over element.
+                average = sum(gradients[v, :, :, elem]/Np^2)
+                gradients[v, :, :, elem] .= average
+            end
+        end
+    end
+
+    return gradients
+end
+
+function compute_flux_array!(Flux, u, direction, equations)
+    Np, nvars = size(Flux)
+    for i in 1:Np
+        Flux[i, :] .= flux(u[i, :], 1, direction, equations)
     end
 end
 
diff --git a/src/equations/compressible_euler_multicomponent_2d.jl b/src/equations/compressible_euler_multicomponent_2d.jl
index 666455fe86c..df088ce304a 100644
--- a/src/equations/compressible_euler_multicomponent_2d.jl
+++ b/src/equations/compressible_euler_multicomponent_2d.jl
@@ -681,6 +681,44 @@ end
     return max(abs(v_ll), abs(v_rr)) + max(c_ll, c_rr)
 end
 
+# Calculate maximum wave speed for local Lax-Friedrichs-type dissipation
+@inline function max_abs_speed_naive(u_ll, u_rr, normal_direction::AbstractVector,
+                                     equations::CompressibleEulerMulticomponentEquations2D)
+    rho_v1_ll, rho_v2_ll, rho_e_ll = u_ll
+    rho_v1_rr, rho_v2_rr, rho_e_rr = u_rr
+
+    # Normalize the normal vector (if not guaranteed to be unit)
+    n1, n2 = normal_direction
+    norm_n = sqrt(n1^2 + n2^2)
+    n1 /= norm_n
+    n2 /= norm_n
+
+    # Densities and gamma
+    rho_ll = density(u_ll, equations)
+    rho_rr = density(u_rr, equations)
+    gamma_ll = totalgamma(u_ll, equations)
+    gamma_rr = totalgamma(u_rr, equations)
+
+    # Velocities
+    v1_ll = rho_v1_ll / rho_ll
+    v2_ll = rho_v2_ll / rho_ll
+    v1_rr = rho_v1_rr / rho_rr
+    v2_rr = rho_v2_rr / rho_rr
+
+    # Normal velocity components
+    v_n_ll = v1_ll * n1 + v2_ll * n2
+    v_n_rr = v1_rr * n1 + v2_rr * n2
+
+    # Pressure and sound speeds
+    p_ll = (gamma_ll - 1f0) * (rho_e_ll - 0.5f0 * rho_ll * (v1_ll^2 + v2_ll^2))
+    p_rr = (gamma_rr - 1f0) * (rho_e_rr - 0.5f0 * rho_rr * (v1_rr^2 + v2_rr^2))
+
+    c_ll = sqrt(gamma_ll * p_ll / rho_ll)
+    c_rr = sqrt(gamma_rr * p_rr / rho_rr)
+
+    return max(abs(v_n_ll), abs(v_n_rr)) + max(c_ll, c_rr)
+end
+
 # Less "cautious", i.e., less overestimating `λ_max` compared to `max_abs_speed_naive`
 @inline function max_abs_speed(u_ll, u_rr, orientation::Integer,
                                equations::CompressibleEulerMulticomponentEquations2D)
@@ -979,4 +1017,22 @@ end
     v = u[orientation] / rho
     return v
 end
+
+@inline function enstrophy_multi_euler(u, gradients,
+                           equations::CompressibleEulerMulticomponentEquations2D)
+    # Enstrophy is 0.5 rho ω⋅ω where ω = ∇ × v
+
+    omega = vorticity(u, gradients, equations)
+
+    return 0.5f0 * omega^2
+end
+
+@inline function vorticity(u, gradients,
+                           equations::CompressibleEulerMulticomponentEquations2D)
+    # Ensure that we have velocity `gradients` by way of the `convert_gradient_variables` function.
+    dv1dx, dv2dx = gradients[1][1:2]
+    dv1dy, dv2dy = gradients[2][1:2]
+
+    return dv2dx - dv1dy
+end
 end # @muladd