diff --git a/src/Trixi.jl b/src/Trixi.jl
index fb307712da8..78ae8ad121f 100644
--- a/src/Trixi.jl
+++ b/src/Trixi.jl
@@ -127,7 +127,8 @@ export TreeMesh
 export DG,
        DGSEM, LobattoLegendreBasis,
        VolumeIntegralWeakForm, VolumeIntegralFluxDifferencing,
-       VolumeIntegralLocalComparison,
+       VolumeIntegralLocalComparison, VolumeIntegralFluxComparison,
+       VolumeIntegralLocalFluxComparison,
        VolumeIntegralPureLGLFiniteVolume,
        VolumeIntegralShockCapturingHG, IndicatorHennemannGassner,
        MortarL2
diff --git a/src/solvers/dg/dg.jl b/src/solvers/dg/dg.jl
index 684a6fb9834..8e0d7b248a1 100644
--- a/src/solvers/dg/dg.jl
+++ b/src/solvers/dg/dg.jl
@@ -91,6 +91,52 @@ function Base.show(io::IO, ::MIME"text/plain", integral::VolumeIntegralLocalComp
 end
 
 
+# TODO: It would be really nice if we could convert other volume integrals to
+#       the flux form required here easily...
+# TODO: This needs a better name...
+struct VolumeIntegralFluxComparison{VolumeFlux_a, VolumeFlux_b} <: AbstractVolumeIntegral
+  volume_flux_a::VolumeFlux_a
+  volume_flux_b::VolumeFlux_b
+end
+
+function Base.show(io::IO, ::MIME"text/plain", integral::VolumeIntegralFluxComparison)
+  @nospecialize integral # reduce precompilation time
+
+  if get(io, :compact, false)
+    show(io, integral)
+  else
+    setup = [
+            "volume flux a" => integral.volume_flux_a,
+            "volume flux b" => integral.volume_flux_b
+            ]
+    summary_box(io, "VolumeIntegralFluxComparison", setup)
+  end
+end
+
+
+# TODO: It would be really nice if we could convert other volume integrals to
+#       the flux form required here easily...
+# TODO: This needs a better name...
+struct VolumeIntegralLocalFluxComparison{VolumeFluxA, VolumeFluxB} <: AbstractVolumeIntegral
+  volume_flux_a::VolumeFluxA
+  volume_flux_b::VolumeFluxB
+end
+
+function Base.show(io::IO, ::MIME"text/plain", integral::VolumeIntegralLocalFluxComparison)
+  @nospecialize integral # reduce precompilation time
+
+  if get(io, :compact, false)
+    show(io, integral)
+  else
+    setup = [
+            "volume flux a" => integral.volume_flux_a,
+            "volume flux b" => integral.volume_flux_b
+            ]
+    summary_box(io, "VolumeIntegralLocalFluxComparison", setup)
+  end
+end
+
+
 """
     VolumeIntegralShockCapturingHG
 
diff --git a/src/solvers/dg/dg_1d.jl b/src/solvers/dg/dg_1d.jl
index dd2040d0791..c6d9879c26a 100644
--- a/src/solvers/dg/dg_1d.jl
+++ b/src/solvers/dg/dg_1d.jl
@@ -60,6 +60,30 @@ function create_cache(mesh::TreeMesh{1}, equations, volume_integral::VolumeInteg
 end
 
 
+function create_cache(mesh::TreeMesh{1}, equations, volume_integral::VolumeIntegralFluxComparison, dg::DG, uEltype)
+
+  have_nonconservative_terms(equations) !== Val(false) && error("Are you kidding me?")
+
+  FluxType = SVector{nvariables(equations), uEltype}
+  w_threaded = [Array{FluxType, 1}(undef, nnodes(dg)) for _ in 1:Threads.nthreads()]
+  fluxes_a_threaded = [Vector{FluxType}(undef, nnodes(dg) - 1) for _ in 1:Threads.nthreads()]
+  fluxes_b_threaded = [Vector{FluxType}(undef, nnodes(dg) - 1) for _ in 1:Threads.nthreads()]
+
+  return (; w_threaded, fluxes_a_threaded, fluxes_b_threaded)
+end
+
+
+function create_cache(mesh::TreeMesh{1}, equations, volume_integral::VolumeIntegralLocalFluxComparison, dg::DG, uEltype)
+
+  have_nonconservative_terms(equations) !== Val(false) && error("Are you kidding me?")
+
+  FluxType = SVector{nvariables(equations), uEltype}
+  w_threaded = [Array{FluxType, 1}(undef, nnodes(dg)) for _ in 1:Threads.nthreads()]
+
+  return (; w_threaded)
+end
+
+
 function create_cache(mesh::TreeMesh{1}, equations,
                       volume_integral::VolumeIntegralShockCapturingHG, dg::DG, uEltype)
   element_ids_dg   = Int[]
@@ -268,6 +292,116 @@ function calc_volume_integral!(du::AbstractArray{<:Any,3}, u,
 end
 
 
+function calc_volume_integral!(du::AbstractArray{<:Any,3}, u,
+                               nonconservative_terms::Val{false}, equations,
+                               volume_integral::VolumeIntegralFluxComparison,
+                               dg::DGSEM, cache)
+  @unpack volume_flux_a, volume_flux_b = volume_integral
+  @unpack inverse_weights = dg.basis
+  @unpack w_threaded, fluxes_a_threaded, fluxes_b_threaded = cache
+
+  @threaded for element in eachelement(dg, cache)
+    w = w_threaded[Threads.threadid()]
+    fluxes_a = fluxes_a_threaded[Threads.threadid()]
+    fluxes_b = fluxes_b_threaded[Threads.threadid()]
+
+    # compute entropy variables
+    for i in eachnode(dg)
+      u_node = get_node_vars(u, equations, dg, i, element)
+      w[i] = cons2entropy(u_node, equations)
+    end
+
+    # compute local high-order fluxes
+    local_fluxes_1!(fluxes_a, u, nonconservative_terms, equations, volume_flux_a, dg, element)
+    local_fluxes_1!(fluxes_b, u, nonconservative_terms, equations, volume_flux_b, dg, element)
+
+    # compare entropy production of both fluxes and choose the more dissipative one
+    for i in eachindex(fluxes_a, fluxes_b)
+      flux_a = fluxes_a[i]
+      flux_b = fluxes_b[i]
+      # if dot(w[i+1] - w[i], flux_a - flux_b) <= 0
+      #   fluxes_a[i] = flux_a # flux_a is more dissipative
+      #   # fluxes_a[i] = 2 * flux_a - flux_b # flux_a is more dissipative
+      # else
+      #   fluxes_a[i] = flux_b # flux_b is more dissipative
+      #   # fluxes_a[i] = 2 * flux_b - flux_a # flux_b is more dissipative
+      # end
+      b = dot(w[i+1] - w[i], flux_b - flux_a)
+      c = 1.0e-7
+      # hyp = sqrt(b^2 + c^2)
+      hyp = hypot(b, c) # sqrt(b^2 + c^2) computed in a numerically stable way
+      # δ = (hyp - b) / hyp # add anti-dissipation as dissipation
+      δ = (hyp - b) / 2hyp # just use the more dissipative flux
+      fluxes_a[i] = flux_a + δ * (flux_b - flux_a)
+    end
+
+    # update volume contribution in locally conservative form
+    add_to_node_vars!(du, inverse_weights[1] * fluxes_a[1], equations, dg, 1, element)
+    for i in 2:nnodes(dg)-1
+      add_to_node_vars!(du, inverse_weights[i] * (fluxes_a[i] - fluxes_a[i-1]), equations, dg, i, element)
+    end
+    add_to_node_vars!(du, -inverse_weights[end] * fluxes_a[end], equations, dg, nnodes(dg), element)
+  end
+end
+
+@inline function local_fluxes_1!(fluxes, u::AbstractArray{<:Any,3},
+                                 nonconservative_terms::Val{false}, equations,
+                                 volume_flux, dg::DGSEM, element)
+  @unpack weights, derivative_split = dg.basis
+
+  for i in eachindex(fluxes)
+    flux1 = zero(eltype(fluxes))
+    for iip in i+1:nnodes(dg)
+      uiip = get_node_vars(u, equations, dg, iip, element)
+      for iim in 1:i
+        uiim = get_node_vars(u, equations, dg, iim, element)
+        flux1 += weights[iim] * derivative_split[iim,iip] * volume_flux(uiim, uiip, 1, equations)
+      end
+    end
+    fluxes[i] = flux1
+  end
+end
+
+
+function calc_volume_integral!(du::AbstractArray{<:Any,3}, u,
+                               nonconservative_terms::Val{false}, equations,
+                               volume_integral::VolumeIntegralLocalFluxComparison,
+                               dg::DGSEM, cache)
+  @unpack volume_flux_a, volume_flux_b = volume_integral
+  @unpack derivative_split = dg.basis
+  @unpack w_threaded = cache
+
+  @threaded for element in eachelement(dg, cache)
+    w = w_threaded[Threads.threadid()]
+
+    # compute entropy variables
+    for i in eachnode(dg)
+      u_node = get_node_vars(u, equations, dg, i, element)
+      w[i] = cons2entropy(u_node, equations)
+    end
+
+    # perform flux-differencing in symmetric form
+    for i in eachnode(dg)
+      u_i = get_node_vars(u, equations, dg, i, element)
+      du_i = zero(u_i)
+      for ii in eachnode(dg)
+        u_ii = get_node_vars(u, equations, dg, ii, element)
+        flux_a = volume_flux_a(u_i, u_ii, 1, equations)
+        flux_b = volume_flux_b(u_i, u_ii, 1, equations)
+        d_i_ii = derivative_split[i, ii]
+        b = -sign(d_i_ii) * dot(w[i] - w[ii], flux_b - flux_a)
+        c = 1.0e-4 # TODO: magic constant determining linear and nonlinear stability 😠
+        hyp = hypot(b, c) # sqrt(b^2 + c^2) computed in a numerically stable way
+        # δ = (hyp - b) / hyp # add anti-dissipation as dissipation
+        δ = (hyp - b) / 2hyp # just use the more dissipative flux
+        du_i += d_i_ii * (flux_a + δ * (flux_b - flux_a))
+      end
+      add_to_node_vars!(du, du_i, equations, dg, i, element)
+    end
+  end
+end
+
+
 # TODO: Taal dimension agnostic
 function calc_volume_integral!(du::AbstractArray{<:Any,3}, u, nonconservative_terms, equations,
                                volume_integral::VolumeIntegralShockCapturingHG,
diff --git a/src/solvers/dg/dg_2d.jl b/src/solvers/dg/dg_2d.jl
index f7da5f8bc2a..9d6b6ce6c76 100644
--- a/src/solvers/dg/dg_2d.jl
+++ b/src/solvers/dg/dg_2d.jl
@@ -75,6 +75,30 @@ function create_cache(mesh::TreeMesh{2}, equations, volume_integral::VolumeInteg
 end
 
 
+function create_cache(mesh::TreeMesh{2}, equations, volume_integral::VolumeIntegralFluxComparison, dg::DG, uEltype)
+
+  have_nonconservative_terms(equations) !== Val(false) && error("Are you kidding me?")
+
+  FluxType = SVector{nvariables(equations), uEltype}
+  w_threaded = [Array{FluxType, 2}(undef, nnodes(dg), nnodes(dg)) for _ in 1:Threads.nthreads()]
+  fluxes_a_threaded = [Vector{FluxType}(undef, nnodes(dg) - 1) for _ in 1:Threads.nthreads()]
+  fluxes_b_threaded = [Vector{FluxType}(undef, nnodes(dg) - 1) for _ in 1:Threads.nthreads()]
+
+  return (; w_threaded, fluxes_a_threaded, fluxes_b_threaded)
+end
+
+
+function create_cache(mesh::TreeMesh{2}, equations, volume_integral::VolumeIntegralLocalFluxComparison, dg::DG, uEltype)
+
+  have_nonconservative_terms(equations) !== Val(false) && error("Are you kidding me?")
+
+  FluxType = SVector{nvariables(equations), uEltype}
+  w_threaded = [Array{FluxType, 2}(undef, nnodes(dg), nnodes(dg)) for _ in 1:Threads.nthreads()]
+
+  return (; w_threaded)
+end
+
+
 function create_cache(mesh::TreeMesh{2}, equations,
                       volume_integral::VolumeIntegralShockCapturingHG, dg::DG, uEltype)
   element_ids_dg   = Int[]
@@ -456,6 +480,172 @@ function calc_volume_integral!(du::AbstractArray{<:Any,4}, u,
 end
 
 
+function calc_volume_integral!(du::AbstractArray{<:Any,4}, u,
+                               nonconservative_terms::Val{false}, equations,
+                               volume_integral::VolumeIntegralFluxComparison,
+                               dg::DGSEM, cache)
+  @unpack volume_flux_a, volume_flux_b = volume_integral
+  @unpack inverse_weights = dg.basis
+  @unpack w_threaded, fluxes_a_threaded, fluxes_b_threaded = cache
+
+  @threaded for element in eachelement(dg, cache)
+    w = w_threaded[Threads.threadid()]
+    fluxes_a = fluxes_a_threaded[Threads.threadid()]
+    fluxes_b = fluxes_b_threaded[Threads.threadid()]
+    c = 1.0e-7
+
+    # compute entropy variables
+    for j in eachnode(dg), i in eachnode(dg)
+      u_node = get_node_vars(u, equations, dg, i, j, element)
+      w[i,j] = cons2entropy(u_node, equations)
+    end
+
+    # x
+    for j in eachnode(dg)
+      # compute local high-order fluxes
+      local_fluxes_1!(fluxes_a, u, nonconservative_terms, equations, volume_flux_a, dg, j, element)
+      local_fluxes_1!(fluxes_b, u, nonconservative_terms, equations, volume_flux_b, dg, j, element)
+
+      # compare entropy production of both fluxes and choose the more dissipative one
+      for i in eachindex(fluxes_a, fluxes_b)
+        flux_a = fluxes_a[i]
+        flux_b = fluxes_b[i]
+        b = dot(w[i+1,j] - w[i,j], flux_b - flux_a)
+        hyp = hypot(b, c) # sqrt(b^2 + c^2) computed in a numerically stable way
+        # δ = (hyp - b) / hyp # add anti-dissipation as dissipation
+        δ = (hyp - b) / 2hyp # just use the more dissipative flux
+        fluxes_a[i] = flux_a + δ * (flux_b - flux_a)
+      end
+
+      # update volume contribution in locally conservative form
+      add_to_node_vars!(du, inverse_weights[1] * fluxes_a[1], equations, dg, 1, j, element)
+      for i in 2:nnodes(dg)-1
+        add_to_node_vars!(du, inverse_weights[i] * (fluxes_a[i] - fluxes_a[i-1]), equations, dg, i, j, element)
+      end
+      add_to_node_vars!(du, -inverse_weights[end] * fluxes_a[end], equations, dg, nnodes(dg), j, element)
+    end
+
+    # y
+    for i in eachnode(dg)
+      # compute local high-order fluxes
+      local_fluxes_2!(fluxes_a, u, nonconservative_terms, equations, volume_flux_a, dg, i, element)
+      local_fluxes_2!(fluxes_b, u, nonconservative_terms, equations, volume_flux_b, dg, i, element)
+
+      # compare entropy production of both fluxes and choose the more dissipative one
+      for j in eachindex(fluxes_a, fluxes_b)
+        flux_a = fluxes_a[j]
+        flux_b = fluxes_b[j]
+        b = dot(w[i,j+1] - w[i,j], flux_b - flux_a)
+        hyp = hypot(b, c) # sqrt(b^2 + c^2) computed in a numerically stable way
+        # δ = (hyp - b) / hyp # add anti-dissipation as dissipation
+        δ = (hyp - b) / 2hyp # just use the more dissipative flux
+        fluxes_a[j] = flux_a + δ * (flux_b - flux_a)
+      end
+
+      # update volume contribution in locally conservative form
+      add_to_node_vars!(du, inverse_weights[1] * fluxes_a[1], equations, dg, i, 1, element)
+      for j in 2:nnodes(dg)-1
+        add_to_node_vars!(du, inverse_weights[j] * (fluxes_a[j] - fluxes_a[j-1]), equations, dg, i, j, element)
+      end
+      add_to_node_vars!(du, -inverse_weights[end] * fluxes_a[end], equations, dg, i, nnodes(dg), element)
+    end
+
+  end
+end
+
+@inline function local_fluxes_1!(fluxes, u::AbstractArray{<:Any,4},
+                                 nonconservative_terms::Val{false}, equations,
+                                 volume_flux, dg::DGSEM, j, element)
+  @unpack weights, derivative_split = dg.basis
+
+  for i in eachindex(fluxes)
+    flux1 = zero(eltype(fluxes))
+    for iip in i+1:nnodes(dg)
+      uiip = get_node_vars(u, equations, dg, iip, j, element)
+      for iim in 1:i
+        uiim = get_node_vars(u, equations, dg, iim, j, element)
+        flux1 += weights[iim] * derivative_split[iim,iip] * volume_flux(uiim, uiip, 1, equations)
+      end
+    end
+    fluxes[i] = flux1
+  end
+end
+
+@inline function local_fluxes_2!(fluxes, u::AbstractArray{<:Any,4},
+                                 nonconservative_terms::Val{false}, equations,
+                                 volume_flux, dg::DGSEM, i, element)
+  @unpack weights, derivative_split = dg.basis
+
+  for j in eachindex(fluxes)
+    flux2 = zero(eltype(fluxes))
+    for jjp in j+1:nnodes(dg)
+      ujjp = get_node_vars(u, equations, dg, i, jjp, element)
+      for jjm in 1:j
+        ujjm = get_node_vars(u, equations, dg, i, jjm, element)
+        flux2 += weights[jjm] * derivative_split[jjm,jjp] * volume_flux(ujjm, ujjp, 2, equations)
+      end
+    end
+    fluxes[j] = flux2
+  end
+end
+
+
+function calc_volume_integral!(du::AbstractArray{<:Any,4}, u,
+                               nonconservative_terms::Val{false}, equations,
+                               volume_integral::VolumeIntegralLocalFluxComparison,
+                               dg::DGSEM, cache)
+  @unpack volume_flux_a, volume_flux_b = volume_integral
+  @unpack derivative_split = dg.basis
+  @unpack w_threaded = cache
+
+  @threaded for element in eachelement(dg, cache)
+    w = w_threaded[Threads.threadid()]
+
+    # compute entropy variables
+    for j in eachnode(dg), i in eachnode(dg)
+      u_node = get_node_vars(u, equations, dg, i, j, element)
+      w[i,j] = cons2entropy(u_node, equations)
+    end
+
+    # perform flux-differencing in symmetric form
+    for j in eachnode(dg), i in eachnode(dg)
+      u_i = get_node_vars(u, equations, dg, i, j, element)
+      du_i = zero(u_i)
+      c = 1.0e-4 # TODO: magic constant determining linear and nonlinear stability 😠
+
+      # x
+      for ii in eachnode(dg)
+        u_ii = get_node_vars(u, equations, dg, ii, j, element)
+        flux_a = volume_flux_a(u_i, u_ii, 1, equations)
+        flux_b = volume_flux_b(u_i, u_ii, 1, equations)
+        d_i_ii = derivative_split[i, ii]
+        b = -sign(d_i_ii) * dot(w[i,j] - w[ii,j], flux_b - flux_a)
+        hyp = hypot(b, c) # sqrt(b^2 + c^2) computed in a numerically stable way
+        # δ = (hyp - b) / hyp # add anti-dissipation as dissipation
+        δ = (hyp - b) / 2hyp # just use the more dissipative flux
+        du_i += d_i_ii * (flux_a + δ * (flux_b - flux_a))
+      end
+
+      # y
+      for jj in eachnode(dg)
+        u_jj = get_node_vars(u, equations, dg, i, jj, element)
+        flux_a = volume_flux_a(u_i, u_jj, 2, equations)
+        flux_b = volume_flux_b(u_i, u_jj, 2, equations)
+        d_j_jj = derivative_split[j, jj]
+        b = -sign(d_j_jj) * dot(w[i,j] - w[i,jj], flux_b - flux_a)
+        hyp = hypot(b, c) # sqrt(b^2 + c^2) computed in a numerically stable way
+        # δ = (hyp - b) / hyp # add anti-dissipation as dissipation
+        δ = (hyp - b) / 2hyp # just use the more dissipative flux
+        du_i += d_j_jj * (flux_a + δ * (flux_b - flux_a))
+      end
+
+      add_to_node_vars!(du, du_i, equations, dg, i, j, element)
+    end
+
+  end
+end
+
+
 # TODO: Taal dimension agnostic
 function calc_volume_integral!(du::AbstractArray{<:Any,4}, u, nonconservative_terms, equations,
                                volume_integral::VolumeIntegralShockCapturingHG,
diff --git a/src/solvers/dg/dg_3d.jl b/src/solvers/dg/dg_3d.jl
index 03c5b6cc973..fb1c6143c79 100644
--- a/src/solvers/dg/dg_3d.jl
+++ b/src/solvers/dg/dg_3d.jl
@@ -94,6 +94,30 @@ function create_cache(mesh::TreeMesh{3}, equations, volume_integral::VolumeInteg
 end
 
 
+function create_cache(mesh::TreeMesh{3}, equations, volume_integral::VolumeIntegralFluxComparison, dg::DG, uEltype)
+
+  have_nonconservative_terms(equations) !== Val(false) && error("Are you kidding me?")
+
+  FluxType = SVector{nvariables(equations), uEltype}
+  w_threaded = [Array{FluxType, 3}(undef, nnodes(dg), nnodes(dg), nnodes(dg)) for _ in 1:Threads.nthreads()]
+  fluxes_a_threaded = [Vector{FluxType}(undef, nnodes(dg) - 1) for _ in 1:Threads.nthreads()]
+  fluxes_b_threaded = [Vector{FluxType}(undef, nnodes(dg) - 1) for _ in 1:Threads.nthreads()]
+
+  return (; w_threaded, fluxes_a_threaded, fluxes_b_threaded)
+end
+
+
+function create_cache(mesh::TreeMesh{3}, equations, volume_integral::VolumeIntegralLocalFluxComparison, dg::DG, uEltype)
+
+  have_nonconservative_terms(equations) !== Val(false) && error("Are you kidding me?")
+
+  FluxType = SVector{nvariables(equations), uEltype}
+  w_threaded = [Array{FluxType, 3}(undef, nnodes(dg), nnodes(dg), nnodes(dg)) for _ in 1:Threads.nthreads()]
+
+  return (; w_threaded)
+end
+
+
 function create_cache(mesh::TreeMesh{3}, equations,
                       volume_integral::VolumeIntegralShockCapturingHG, dg::DG, uEltype)
   element_ids_dg   = Int[]
@@ -476,6 +500,228 @@ function calc_volume_integral!(du::AbstractArray{<:Any,5}, u,
 end
 
 
+function calc_volume_integral!(du::AbstractArray{<:Any,5}, u,
+                               nonconservative_terms::Val{false}, equations,
+                               volume_integral::VolumeIntegralFluxComparison,
+                               dg::DGSEM, cache)
+  @unpack volume_flux_a, volume_flux_b = volume_integral
+  @unpack inverse_weights = dg.basis
+  @unpack w_threaded, fluxes_a_threaded, fluxes_b_threaded = cache
+
+  @threaded for element in eachelement(dg, cache)
+    w = w_threaded[Threads.threadid()]
+    fluxes_a = fluxes_a_threaded[Threads.threadid()]
+    fluxes_b = fluxes_b_threaded[Threads.threadid()]
+    c = 1.0e-7
+
+    # compute entropy variables
+    for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg)
+      u_node = get_node_vars(u, equations, dg, i, j, k, element)
+      w[i,j,k] = cons2entropy(u_node, equations)
+    end
+
+    # x
+    for k in eachnode(dg), j in eachnode(dg)
+      # compute local high-order fluxes
+      local_fluxes_1!(fluxes_a, u, nonconservative_terms, equations, volume_flux_a, dg, j, k, element)
+      local_fluxes_1!(fluxes_b, u, nonconservative_terms, equations, volume_flux_b, dg, j, k, element)
+
+      # compare entropy production of both fluxes and choose the more dissipative one
+      for i in eachindex(fluxes_a, fluxes_b)
+        flux_a = fluxes_a[i]
+        flux_b = fluxes_b[i]
+        b = dot(w[i+1,j,k] - w[i,j,k], flux_b - flux_a)
+        hyp = hypot(b, c) # sqrt(b^2 + c^2) computed in a numerically stable way
+        # δ = (hyp - b) / hyp # add anti-dissipation as dissipation
+        δ = (hyp - b) / 2hyp # just use the more dissipative flux
+        fluxes_a[i] = flux_a + δ * (flux_b - flux_a)
+      end
+
+      # update volume contribution in locally conservative form
+      add_to_node_vars!(du, inverse_weights[1] * fluxes_a[1], equations, dg, 1, j, k, element)
+      for i in 2:nnodes(dg)-1
+        add_to_node_vars!(du, inverse_weights[i] * (fluxes_a[i] - fluxes_a[i-1]), equations, dg, i, j, k, element)
+      end
+      add_to_node_vars!(du, -inverse_weights[end] * fluxes_a[end], equations, dg, nnodes(dg), j, k, element)
+    end
+
+    # y
+    for k in eachnode(dg), i in eachnode(dg)
+      # compute local high-order fluxes
+      local_fluxes_2!(fluxes_a, u, nonconservative_terms, equations, volume_flux_a, dg, i, k, element)
+      local_fluxes_2!(fluxes_b, u, nonconservative_terms, equations, volume_flux_b, dg, i, k, element)
+
+      # compare entropy production of both fluxes and choose the more dissipative one
+      for j in eachindex(fluxes_a, fluxes_b)
+        flux_a = fluxes_a[j]
+        flux_b = fluxes_b[j]
+        b = dot(w[i,j+1,k] - w[i,j,k], flux_b - flux_a)
+        hyp = hypot(b, c) # sqrt(b^2 + c^2) computed in a numerically stable way
+        # δ = (hyp - b) / hyp # add anti-dissipation as dissipation
+        δ = (hyp - b) / 2hyp # just use the more dissipative flux
+        fluxes_a[j] = flux_a + δ * (flux_b - flux_a)
+      end
+
+      # update volume contribution in locally conservative form
+      add_to_node_vars!(du, inverse_weights[1] * fluxes_a[1], equations, dg, i, 1, k, element)
+      for j in 2:nnodes(dg)-1
+        add_to_node_vars!(du, inverse_weights[j] * (fluxes_a[j] - fluxes_a[j-1]), equations, dg, i, j, k, element)
+      end
+      add_to_node_vars!(du, -inverse_weights[end] * fluxes_a[end], equations, dg, i, nnodes(dg), k, element)
+    end
+
+    # z
+    for j in eachnode(dg), i in eachnode(dg)
+      # compute local high-order fluxes
+      local_fluxes_3!(fluxes_a, u, nonconservative_terms, equations, volume_flux_a, dg, i, j, element)
+      local_fluxes_3!(fluxes_b, u, nonconservative_terms, equations, volume_flux_b, dg, i, j, element)
+
+      # compare entropy production of both fluxes and choose the more dissipative one
+      for k in eachindex(fluxes_a, fluxes_b)
+        flux_a = fluxes_a[k]
+        flux_b = fluxes_b[k]
+        b = dot(w[i,j,k+1] - w[i,j,k], flux_b - flux_a)
+        hyp = hypot(b, c) # sqrt(b^2 + c^2) computed in a numerically stable way
+        # δ = (hyp - b) / hyp # add anti-dissipation as dissipation
+        δ = (hyp - b) / 2hyp # just use the more dissipative flux
+        fluxes_a[k] = flux_a + δ * (flux_b - flux_a)
+      end
+
+      # update volume contribution in locally conservative form
+      add_to_node_vars!(du, inverse_weights[1] * fluxes_a[1], equations, dg, i, j, 1, element)
+      for k in 2:nnodes(dg)-1
+        add_to_node_vars!(du, inverse_weights[k] * (fluxes_a[k] - fluxes_a[k-1]), equations, dg, i, j, k, element)
+      end
+      add_to_node_vars!(du, -inverse_weights[end] * fluxes_a[end], equations, dg, i, j, nnodes(dg), element)
+    end
+
+  end
+end
+
+@inline function local_fluxes_1!(fluxes, u::AbstractArray{<:Any,5},
+                                 nonconservative_terms::Val{false}, equations,
+                                 volume_flux, dg::DGSEM, j, k, element)
+  @unpack weights, derivative_split = dg.basis
+
+  for i in eachindex(fluxes)
+    flux1 = zero(eltype(fluxes))
+    for iip in i+1:nnodes(dg)
+      uiip = get_node_vars(u, equations, dg, iip, j, k, element)
+      for iim in 1:i
+        uiim = get_node_vars(u, equations, dg, iim, j, k, element)
+        flux1 += weights[iim] * derivative_split[iim,iip] * volume_flux(uiim, uiip, 1, equations)
+      end
+    end
+    fluxes[i] = flux1
+  end
+end
+
+@inline function local_fluxes_2!(fluxes, u::AbstractArray{<:Any,5},
+                                 nonconservative_terms::Val{false}, equations,
+                                 volume_flux, dg::DGSEM, i, k, element)
+  @unpack weights, derivative_split = dg.basis
+
+  for j in eachindex(fluxes)
+    flux2 = zero(eltype(fluxes))
+    for jjp in j+1:nnodes(dg)
+      ujjp = get_node_vars(u, equations, dg, i, jjp, k, element)
+      for jjm in 1:j
+        ujjm = get_node_vars(u, equations, dg, i, jjm, k, element)
+        flux2 += weights[jjm] * derivative_split[jjm,jjp] * volume_flux(ujjm, ujjp, 2, equations)
+      end
+    end
+    fluxes[j] = flux2
+  end
+end
+
+@inline function local_fluxes_3!(fluxes, u::AbstractArray{<:Any,5},
+                                 nonconservative_terms::Val{false}, equations,
+                                 volume_flux, dg::DGSEM, i, j, element)
+  @unpack weights, derivative_split = dg.basis
+
+  for k in eachindex(fluxes)
+    flux3 = zero(eltype(fluxes))
+    for kkp in k+1:nnodes(dg)
+      ukkp = get_node_vars(u, equations, dg, i, j, kkp, element)
+      for kkm in 1:k
+        ukkm = get_node_vars(u, equations, dg, i, j, kkm, element)
+        flux3 += weights[kkm] * derivative_split[kkm,kkp] * volume_flux(ukkm, ukkp, 3, equations)
+      end
+    end
+    fluxes[k] = flux3
+  end
+end
+
+
+function calc_volume_integral!(du::AbstractArray{<:Any,5}, u,
+                               nonconservative_terms::Val{false}, equations,
+                               volume_integral::VolumeIntegralLocalFluxComparison,
+                               dg::DGSEM, cache)
+  @unpack volume_flux_a, volume_flux_b = volume_integral
+  @unpack derivative_split = dg.basis
+  @unpack w_threaded = cache
+
+  @threaded for element in eachelement(dg, cache)
+    w = w_threaded[Threads.threadid()]
+
+    # compute entropy variables
+    for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg)
+      u_node = get_node_vars(u, equations, dg, i, j, k, element)
+      w[i,j,k] = cons2entropy(u_node, equations)
+    end
+
+    # perform flux-differencing in symmetric form
+    for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg)
+      u_i = get_node_vars(u, equations, dg, i, j, k, element)
+      du_i = zero(u_i)
+      c = 1.0e-4 # TODO: magic constant determining linear and nonlinear stability 😠
+
+      # x
+      for ii in eachnode(dg)
+        u_ii = get_node_vars(u, equations, dg, ii, j, k, element)
+        flux_a = volume_flux_a(u_i, u_ii, 1, equations)
+        flux_b = volume_flux_b(u_i, u_ii, 1, equations)
+        d_i_ii = derivative_split[i, ii]
+        b = -sign(d_i_ii) * dot(w[i,j,k] - w[ii,j,k], flux_b - flux_a)
+        hyp = hypot(b, c) # sqrt(b^2 + c^2) computed in a numerically stable way
+        # δ = (hyp - b) / hyp # add anti-dissipation as dissipation
+        δ = (hyp - b) / 2hyp # just use the more dissipative flux
+        du_i += d_i_ii * (flux_a + δ * (flux_b - flux_a))
+      end
+
+      # y
+      for jj in eachnode(dg)
+        u_jj = get_node_vars(u, equations, dg, i, jj, k, element)
+        flux_a = volume_flux_a(u_i, u_jj, 2, equations)
+        flux_b = volume_flux_b(u_i, u_jj, 2, equations)
+        d_j_jj = derivative_split[j, jj]
+        b = -sign(d_j_jj) * dot(w[i,j,k] - w[i,jj,k], flux_b - flux_a)
+        hyp = hypot(b, c) # sqrt(b^2 + c^2) computed in a numerically stable way
+        # δ = (hyp - b) / hyp # add anti-dissipation as dissipation
+        δ = (hyp - b) / 2hyp # just use the more dissipative flux
+        du_i += d_j_jj * (flux_a + δ * (flux_b - flux_a))
+      end
+
+      # z
+      for kk in eachnode(dg)
+        u_kk = get_node_vars(u, equations, dg, i, j, kk, element)
+        flux_a = volume_flux_a(u_i, u_kk, 3, equations)
+        flux_b = volume_flux_b(u_i, u_kk, 3, equations)
+        d_k_kk = derivative_split[k, kk]
+        b = -sign(d_k_kk) * dot(w[i,j,k] - w[i,j,kk], flux_b - flux_a)
+        hyp = hypot(b, c) # sqrt(b^2 + c^2) computed in a numerically stable way
+        # δ = (hyp - b) / hyp # add anti-dissipation as dissipation
+        δ = (hyp - b) / 2hyp # just use the more dissipative flux
+        du_i += d_k_kk * (flux_a + δ * (flux_b - flux_a))
+      end
+
+      add_to_node_vars!(du, du_i, equations, dg, i, j, k, element)
+    end
+
+  end
+end
+
+
 # TODO: Taal dimension agnostic
 function calc_volume_integral!(du::AbstractArray{<:Any,5}, u, nonconservative_terms, equations,
                                volume_integral::VolumeIntegralShockCapturingHG,