diff --git a/examples/2d/elixir_euler_kelvin_helmholtz_instability_comparison_pure_ec.jl b/examples/2d/elixir_euler_kelvin_helmholtz_instability_comparison_pure_ec.jl
new file mode 100644
index 00000000000..71b0cdfc366
--- /dev/null
+++ b/examples/2d/elixir_euler_kelvin_helmholtz_instability_comparison_pure_ec.jl
@@ -0,0 +1,62 @@
+
+# TODO: Needs tests
+using Random: seed!
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the compressible Euler equations
+gamma = 1.4
+equations = CompressibleEulerEquations2D(gamma)
+
+seed!(0)
+initial_condition = initial_condition_khi
+
+surface_flux = flux_lax_friedrichs
+volume_flux  = flux_chandrashekar
+polydeg = 3
+basis = LobattoLegendreBasis(polydeg)
+volume_integral = VolumeIntegralLocalComparison(VolumeIntegralFluxDifferencing(volume_flux),
+                                                (;alpha_ec=1, alpha_cen=0), variant=:parametric)
+solver = DGSEM(3, surface_flux, volume_integral)
+
+coordinates_min = (-0.5, -0.5)
+coordinates_max = ( 0.5,  0.5)
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level=5,
+                n_cells_max=100_000)
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 5.0)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 100
+analysis_callback = AnalysisCallback(semi, interval=analysis_interval)
+
+alive_callback = AliveCallback(analysis_interval=analysis_interval)
+
+save_solution = SaveSolutionCallback(interval=20,
+                                     save_initial_solution=true,
+                                     save_final_solution=true,
+                                     solution_variables=cons2prim)
+
+stepsize_callback = StepsizeCallback(cfl=1.3)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback, alive_callback,
+                        save_solution,
+                        stepsize_callback)
+
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
+            dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
+            save_everystep=false, callback=callbacks);
+summary_callback() # print the timer summary
diff --git a/src/solvers/dg/dg.jl b/src/solvers/dg/dg.jl
index b01db513dc7..684a6fb9834 100644
--- a/src/solvers/dg/dg.jl
+++ b/src/solvers/dg/dg.jl
@@ -63,18 +63,28 @@ end
 #       The choice coded below would then be the special case of blending the
 #        weak form with flux differencing, but we could also use other choices.
 # TODO: This needs a better name...
-struct VolumeIntegralLocalComparison{VolumeFlux} <: AbstractVolumeIntegral
+struct VolumeIntegralLocalComparison{Variant, VolumeFlux, Parameters} <: AbstractVolumeIntegral
   volume_integral_flux_differencing::VolumeIntegralFluxDifferencing{VolumeFlux}
+  parameters::Parameters
 end
 
-function Base.show(io::IO, ::MIME"text/plain", integral::VolumeIntegralLocalComparison)
+function VolumeIntegralLocalComparison(volume_integral_flux_differencing, parameters=NamedTuple();
+                                       variant=:default)
+  VolumeIntegralLocalComparison{Val{variant},
+                                typeof(volume_integral_flux_differencing.volume_flux),
+                                typeof(parameters)}(volume_integral_flux_differencing, parameters)
+end
+
+function Base.show(io::IO, ::MIME"text/plain", integral::VolumeIntegralLocalComparison{Variant}) where Variant
   @nospecialize integral # reduce precompilation time
 
   if get(io, :compact, false)
     show(io, integral)
   else
     setup = [
-            "volume flux" => integral.volume_integral_flux_differencing.volume_flux
+            "variant" => Variant,
+            "volume flux" => integral.volume_integral_flux_differencing.volume_flux,
+            "parameters" => integral.parameters
             ]
     summary_box(io, "VolumeIntegralLocalComparison", setup)
   end
diff --git a/src/solvers/dg/dg_2d.jl b/src/solvers/dg/dg_2d.jl
index c1dbc03ab1e..f7da5f8bc2a 100644
--- a/src/solvers/dg/dg_2d.jl
+++ b/src/solvers/dg/dg_2d.jl
@@ -374,7 +374,7 @@ end
 
 function calc_volume_integral!(du::AbstractArray{<:Any,4}, u,
                                nonconservative_terms, equations,
-                               volume_integral::VolumeIntegralLocalComparison,
+                               volume_integral::VolumeIntegralLocalComparison{Val{:default}},
                                dg::DGSEM, cache)
   @unpack weights = dg.basis
   @unpack du_ec, du_cen = cache
@@ -412,6 +412,50 @@ function calc_volume_integral!(du::AbstractArray{<:Any,4}, u,
 end
 
 
+function calc_volume_integral!(du::AbstractArray{<:Any,4}, u,
+                               nonconservative_terms, equations,
+                               volume_integral::VolumeIntegralLocalComparison{Val{:parametric}},
+                               dg::DGSEM, cache)
+  @unpack weights = dg.basis
+  @unpack du_ec, du_cen = cache
+  @unpack volume_flux = volume_integral.volume_integral_flux_differencing
+  @unpack parameters = volume_integral
+  @unpack alpha_ec = parameters
+  @unpack alpha_cen = parameters
+
+  du_ec  .= zero(eltype(du_ec))
+  du_cen .= zero(eltype(du_cen))
+
+  @threaded for element in eachelement(dg, cache)
+    # compute volume integral with flux, and for comparison with central flux
+    split_form_kernel!(du_ec,  u, nonconservative_terms, equations, volume_flux, dg, cache, element)
+    weak_form_kernel!(du_cen, u, nonconservative_terms, equations, dg, cache, element)
+
+    # compute entropy production of both volume integrals
+    delta_entropy = zero(eltype(du))
+    for j in eachnode(dg), i in eachnode(dg)
+      du_ec_node  = get_node_vars(du_ec,  equations, dg, i, j, element)
+      du_cen_node = get_node_vars(du_cen, equations, dg, i, j, element)
+      w_node = cons2entropy(get_node_vars(u, equations, dg, i, j, element), equations)
+      delta_entropy += weights[i] * weights[j] * dot(w_node, du_ec_node - du_cen_node)
+    end
+    if delta_entropy < 0
+      for j in eachnode(dg), i in eachnode(dg)
+        du_cen_node = get_node_vars(du_cen, equations, dg, i, j, element)
+        add_to_node_vars!(du, du_cen_node, equations, dg, i, j, element)
+      end
+    else
+      for j in eachnode(dg), i in eachnode(dg)
+        du_ec_node = get_node_vars(du_ec, equations, dg, i, j, element)
+        du_cen_node = get_node_vars(du_cen, equations, dg, i, j, element)
+        add_to_node_vars!(du, alpha_ec*du_ec_node + alpha_cen * du_cen_node,
+                          equations, dg, i, j, element)
+      end
+    end
+  end
+end
+
+
 # TODO: Taal dimension agnostic
 function calc_volume_integral!(du::AbstractArray{<:Any,4}, u, nonconservative_terms, equations,
                                volume_integral::VolumeIntegralShockCapturingHG,