diff --git a/examples/structured_1d_dgsem/elixir_traffic_flow_lwr_greenlight.jl b/examples/structured_1d_dgsem/elixir_traffic_flow_lwr_greenlight.jl
index ea3102ab1e2..34cc7f0fba0 100644
--- a/examples/structured_1d_dgsem/elixir_traffic_flow_lwr_greenlight.jl
+++ b/examples/structured_1d_dgsem/elixir_traffic_flow_lwr_greenlight.jl
@@ -7,9 +7,7 @@ using Trixi
 equations = TrafficFlowLWREquations1D()
 
 basis = LobattoLegendreBasis(3)
-
 surface_flux = flux_hll
-                                                 
 solver = DGSEM(basis, surface_flux)
 
 coordinates_min = (-1.0,) # minimum coordinate
@@ -17,15 +15,16 @@ coordinates_max = (1.0,) # maximum coordinate
 cells_per_dimension = (64,)
 
 # Create curved mesh with 16 cells
-mesh = StructuredMesh(cells_per_dimension, coordinates_min, coordinates_max, periodicity = false)
+mesh = StructuredMesh(cells_per_dimension, coordinates_min, coordinates_max,
+                      periodicity = false)
 
 # Example inspired from http://www.clawpack.org/riemann_book/html/Traffic_flow.html#Example:-green-light
 # Green light that at x = 0 which switches at t = 0 from red to green.
 # To the left there are cars bumper to bumper, to the right there are no cars.
 function initial_condition_greenlight(x, t, equation::TrafficFlowLWREquations1D)
-  scalar = x[1] < 0.0 ? 1.0 : 0.0
+    scalar = x[1] < 0.0 ? 1.0 : 0.0
 
-  return SVector(scalar)
+    return SVector(scalar)
 end
 
 ###############################################################################
@@ -33,28 +32,27 @@ end
 
 # Assume that there are always cars waiting at the left 
 function inflow(x, t, equations::TrafficFlowLWREquations1D)
-  return initial_condition_greenlight(coordinate_min, t, equations)
+    return initial_condition_greenlight(coordinate_min, t, equations)
 end
 boundary_condition_inflow = BoundaryConditionDirichlet(inflow)
 
 # Cars may leave the modeled domain
 function boundary_condition_outflow(u_inner, orientation, normal_direction, x, t,
-                                   surface_flux_function, equations::TrafficFlowLWREquations1D)
-  # Calculate the boundary flux entirely from the internal solution state
-  flux = Trixi.flux(u_inner, orientation, equations)
+                                    surface_flux_function,
+                                    equations::TrafficFlowLWREquations1D)
+    # Calculate the boundary flux entirely from the internal solution state
+    flux = Trixi.flux(u_inner, orientation, equations)
 
-  return flux
+    return flux
 end
 
+boundary_conditions = (x_neg = boundary_condition_outflow,
+                       x_pos = boundary_condition_outflow)
 
-boundary_conditions = (x_neg=boundary_condition_outflow,
-                       x_pos=boundary_condition_outflow)
-                       
 initial_condition = initial_condition_greenlight
 
 semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
-                                    boundary_conditions=boundary_conditions)
-
+                                    boundary_conditions = boundary_conditions)
 
 ###############################################################################
 # ODE solvers, callbacks etc.
@@ -65,12 +63,11 @@ ode = semidiscretize(semi, tspan)
 summary_callback = SummaryCallback()
 
 analysis_interval = 100
-analysis_callback = AnalysisCallback(semi, interval=analysis_interval)
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
 
-alive_callback = AliveCallback(analysis_interval=analysis_interval)
-
-stepsize_callback = StepsizeCallback(cfl=1.2)
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
 
+stepsize_callback = StepsizeCallback(cfl = 1.2)
 
 callbacks = CallbackSet(summary_callback,
                         analysis_callback, alive_callback,
@@ -79,9 +76,8 @@ callbacks = CallbackSet(summary_callback,
 ###############################################################################
 # run the simulation
 
-
-sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
-            dt=42, # solve needs some value here but it will be overwritten by the stepsize_callback
-            save_everystep=false, callback=callbacks);
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false),
+            dt = 42, # 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/examples/tree_1d_dgsem/elixir_traffic_flow_lwr_convergence.jl b/examples/tree_1d_dgsem/elixir_traffic_flow_lwr_convergence.jl
index 3153f9977ab..f887a198023 100644
--- a/examples/tree_1d_dgsem/elixir_traffic_flow_lwr_convergence.jl
+++ b/examples/tree_1d_dgsem/elixir_traffic_flow_lwr_convergence.jl
@@ -5,7 +5,7 @@ using Trixi
 ###############################################################################
 
 equations = TrafficFlowLWREquations1D()
-   
+
 # Use first order finite volume to prevent oscillations at the shock
 solver = DGSEM(polydeg = 3, surface_flux = flux_hll)
 
@@ -24,7 +24,6 @@ initial_condition = initial_condition_convergence_test
 semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
                                     source_terms = source_terms_convergence_test)
 
-
 ###############################################################################
 # ODE solvers, callbacks etc.
 
@@ -34,24 +33,22 @@ ode = semidiscretize(semi, tspan)
 summary_callback = SummaryCallback()
 
 analysis_interval = 100
-analysis_callback = AnalysisCallback(semi, interval=analysis_interval)
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
 
-alive_callback = AliveCallback(analysis_interval=analysis_interval)
-
-stepsize_callback = StepsizeCallback(cfl=1.6)
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
 
+stepsize_callback = StepsizeCallback(cfl = 1.6)
 
 callbacks = CallbackSet(summary_callback,
-                        analysis_callback, 
+                        analysis_callback,
                         alive_callback,
                         stepsize_callback)
 
 ###############################################################################
 # run the simulation
 
-
-sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
-            dt=42, # solve needs some value here but it will be overwritten by the stepsize_callback
-            save_everystep=false, callback=callbacks);
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false),
+            dt = 42, # 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/examples/tree_1d_dgsem/elixir_traffic_flow_lwr_trafficjam.jl b/examples/tree_1d_dgsem/elixir_traffic_flow_lwr_trafficjam.jl
index b2529026930..d92d91dc1b6 100644
--- a/examples/tree_1d_dgsem/elixir_traffic_flow_lwr_trafficjam.jl
+++ b/examples/tree_1d_dgsem/elixir_traffic_flow_lwr_trafficjam.jl
@@ -5,7 +5,7 @@ using Trixi
 ###############################################################################
 
 equations = TrafficFlowLWREquations1D()
-   
+
 # Use first order finite volume to prevent oscillations at the shock
 solver = DGSEM(polydeg = 0, surface_flux = flux_lax_friedrichs)
 
@@ -22,35 +22,35 @@ mesh = TreeMesh(coordinates_min, coordinates_max,
 # Discontinuous initial condition (Riemann Problem) leading to a shock that moves to the left.
 # The shock corresponds to the traffic congestion.
 function initial_condition_traffic_jam(x, t, equation::TrafficFlowLWREquations1D)
-  scalar = x[1] < 0.0 ? 0.5 : 1.0
+    scalar = x[1] < 0.0 ? 0.5 : 1.0
 
-  return SVector(scalar)
+    return SVector(scalar)
 end
 
 ###############################################################################
 # Specify non-periodic boundary conditions
 
 function outflow(x, t, equations::TrafficFlowLWREquations1D)
-  return initial_condition_traffic_jam(coordinates_min, t, equations)
+    return initial_condition_traffic_jam(coordinates_min, t, equations)
 end
 boundary_condition_outflow = BoundaryConditionDirichlet(outflow)
 
 function boundary_condition_inflow(u_inner, orientation, normal_direction, x, t,
-                                   surface_flux_function, equations::TrafficFlowLWREquations1D)
-  # Calculate the boundary flux entirely from the internal solution state
-  flux = Trixi.flux(u_inner, orientation, equations)
+                                   surface_flux_function,
+                                   equations::TrafficFlowLWREquations1D)
+    # Calculate the boundary flux entirely from the internal solution state
+    flux = Trixi.flux(u_inner, orientation, equations)
 
-  return flux
+    return flux
 end
 
-boundary_conditions = (x_neg=boundary_condition_outflow,
-                       x_pos=boundary_condition_inflow)
-                       
+boundary_conditions = (x_neg = boundary_condition_outflow,
+                       x_pos = boundary_condition_inflow)
+
 initial_condition = initial_condition_traffic_jam
 
 semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
-                                    boundary_conditions=boundary_conditions)
-
+                                    boundary_conditions = boundary_conditions)
 
 ###############################################################################
 # ODE solvers, callbacks etc.
@@ -61,12 +61,11 @@ ode = semidiscretize(semi, tspan)
 summary_callback = SummaryCallback()
 
 analysis_interval = 100
-analysis_callback = AnalysisCallback(semi, interval=analysis_interval)
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
 
-alive_callback = AliveCallback(analysis_interval=analysis_interval)
-
-stepsize_callback = StepsizeCallback(cfl=2.0)
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
 
+stepsize_callback = StepsizeCallback(cfl = 2.0)
 
 callbacks = CallbackSet(summary_callback,
                         analysis_callback, alive_callback,
@@ -75,9 +74,8 @@ callbacks = CallbackSet(summary_callback,
 ###############################################################################
 # run the simulation
 
-
-sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
-            dt=42, # solve needs some value here but it will be overwritten by the stepsize_callback
-            save_everystep=false, callback=callbacks);
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false),
+            dt = 42, # 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/equations/equations.jl b/src/equations/equations.jl
index 62f2f1d2c27..65875a2a7e5 100644
--- a/src/equations/equations.jl
+++ b/src/equations/equations.jl
@@ -509,7 +509,7 @@ abstract type AbstractEquationsParabolic{NDIMS, NVARS, GradientVariables} <:
               AbstractEquations{NDIMS, NVARS} end
 
 # Lighthill-Witham-Richards (LWR) traffic flow model
-abstract type AbstractTrafficFlowLWREquations{NDIMS, NVARS} <: 
+abstract type AbstractTrafficFlowLWREquations{NDIMS, NVARS} <:
               AbstractEquations{NDIMS, NVARS} end
-include("traffic_flow_lwr_1d.jl")              
+include("traffic_flow_lwr_1d.jl")
 end # @muladd
diff --git a/src/equations/traffic_flow_lwr_1d.jl b/src/equations/traffic_flow_lwr_1d.jl
index 59d82f68700..f2d4836ce90 100644
--- a/src/equations/traffic_flow_lwr_1d.jl
+++ b/src/equations/traffic_flow_lwr_1d.jl
@@ -3,120 +3,114 @@
 # we need to opt-in explicitly.
 # See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
 @muladd begin
-
-
-@doc raw"""
-    TrafficFlowLWREquations1D
-
-The classic Lighthill-Witham Richards (LWR) model for 1D traffic flow.
-```math
-\partial_t u + v_{\text{max}} \partial_1 [u (1 - u)] = 0
-```
-See e.g. Section 11.1 of 
-- Randall LeVeque (2002)
-Finite Volume Methods for Hyperbolic Problems
-[DOI: 10.1017/CBO9780511791253]https://doi.org/10.1017/CBO9780511791253
-"""
-struct TrafficFlowLWREquations1D{RealT<:Real} <: AbstractTrafficFlowLWREquations{1, 1} 
-  v_max::RealT
-
-  function TrafficFlowLWREquations1D(v_max=1.0)
-    new{typeof(v_max)}(v_max)
-  end
-end
-
-
-varnames(::typeof(cons2cons), ::TrafficFlowLWREquations1D) = ("car-density", )
-varnames(::typeof(cons2prim), ::TrafficFlowLWREquations1D) = ("car-density", )
-
-
-"""
-    initial_condition_constant(x, t, equations::TrafficFlowLWREquations1D)
-
-A constant initial condition to test free-stream preservation.
-"""
-function initial_condition_constant(x, t, equations::TrafficFlowLWREquations1D)
-  return SVector(1.0)
-end
-
-
-"""
-    initial_condition_convergence_test(x, t, equations::TrafficFlowLWREquations1D)
-
-A smooth initial condition used for convergence tests.
-"""
-function initial_condition_convergence_test(x, t, equations::TrafficFlowLWREquations1D)
-  c = 2.0
-  A = 1.0
-  L = 1
-  f = 1/L
-  omega = 2 * pi * f
-  scalar = c + A * sin(omega * (x[1] - t))
-
-  return SVector(scalar)
-end
-
-
-"""
-    source_terms_convergence_test(u, x, t, equations::TrafficFlowLWREquations1D)
-
-Source terms used for convergence tests in combination with
-[`initial_condition_convergence_test`](@ref).
-"""
-@inline function source_terms_convergence_test(u, x, t, equations::TrafficFlowLWREquations1D)
-  # Same settings as in `initial_condition`
-  c = 2.0
-  A = 1.0
-  L = 1
-  f = 1/L
-  omega = 2 * pi * f
-  du = omega * cos(omega * (x[1] - t)) * (-1 -2 * equations.v_max * sin(omega * (x[1] - t)) - 3 * equations.v_max)
-
-  return SVector(du)
-end
-
-
-# Calculate 1D flux in for a single point
-@inline function flux(u, orientation::Integer, equations::TrafficFlowLWREquations1D)
-  return SVector(equations.v_max * u[1] * (1.0 - u[1]))
-end
-
-# Calculate maximum wave speed for local Lax-Friedrichs-type dissipation
-@inline function max_abs_speed_naive(u_ll, u_rr, orientation::Integer, equations::TrafficFlowLWREquations1D)
-  λ_max = max(abs(equations.v_max * (1.0 - 2 * u_ll[1])), 
-              abs(equations.v_max * (1.0 - 2 * u_rr[1])))
-end
-
-# Calculate minimum and maximum wave speeds for HLL-type fluxes
-@inline function min_max_speed_naive(u_ll, u_rr, orientation::Integer, equations::TrafficFlowLWREquations1D)
-  jac_L = equations.v_max * (1.0 - 2 * u_ll[1])
-  jac_R = equations.v_max * (1.0 - 2 * u_rr[1])
-  
-  λ_min = min(jac_L, jac_R)
-  λ_max = max(jac_L, jac_R)
-
-  return λ_min, λ_max
-end
-
-@inline function max_abs_speeds(u, equations::TrafficFlowLWREquations1D)
-  return (abs(equations.v_max * (1.0 - 2 * u[1])),)
-end
-
-# Convert conservative variables to primitive
-@inline cons2prim(u, equations::TrafficFlowLWREquations1D) = u
-
-# Convert conservative variables to entropy variables
-@inline cons2entropy(u, equations::TrafficFlowLWREquations1D) = u
-
-
-# Calculate entropy for a conservative state `cons`
-@inline entropy(u::Real, ::TrafficFlowLWREquations1D) = 0.5 * u^2
-@inline entropy(u, equations::TrafficFlowLWREquations1D) = entropy(u[1], equations)
-
-
-# Calculate total energy for a conservative state `cons`
-@inline energy_total(u::Real, ::TrafficFlowLWREquations1D) = 0.5 * u^2
-@inline energy_total(u, equations::TrafficFlowLWREquations1D) = energy_total(u[1], equations)
-
-
+    @doc raw"""
+        TrafficFlowLWREquations1D
+
+    The classic Lighthill-Witham Richards (LWR) model for 1D traffic flow.
+    ```math
+    \partial_t u + v_{\text{max}} \partial_1 [u (1 - u)] = 0
+    ```
+    See e.g. Section 11.1 of 
+    - Randall LeVeque (2002)
+    Finite Volume Methods for Hyperbolic Problems
+    [DOI: 10.1017/CBO9780511791253]https://doi.org/10.1017/CBO9780511791253
+    """
+    struct TrafficFlowLWREquations1D{RealT <: Real} <: AbstractTrafficFlowLWREquations{1, 1}
+        v_max::RealT
+
+        function TrafficFlowLWREquations1D(v_max = 1.0)
+            new{typeof(v_max)}(v_max)
+        end
+    end
+
+    varnames(::typeof(cons2cons), ::TrafficFlowLWREquations1D) = ("car-density",)
+    varnames(::typeof(cons2prim), ::TrafficFlowLWREquations1D) = ("car-density",)
+
+    """
+        initial_condition_constant(x, t, equations::TrafficFlowLWREquations1D)
+
+    A constant initial condition to test free-stream preservation.
+    """
+    function initial_condition_constant(x, t, equations::TrafficFlowLWREquations1D)
+        return SVector(1.0)
+    end
+
+    """
+        initial_condition_convergence_test(x, t, equations::TrafficFlowLWREquations1D)
+
+    A smooth initial condition used for convergence tests.
+    """
+    function initial_condition_convergence_test(x, t, equations::TrafficFlowLWREquations1D)
+        c = 2.0
+        A = 1.0
+        L = 1
+        f = 1 / L
+        omega = 2 * pi * f
+        scalar = c + A * sin(omega * (x[1] - t))
+
+        return SVector(scalar)
+    end
+
+    """
+        source_terms_convergence_test(u, x, t, equations::TrafficFlowLWREquations1D)
+
+    Source terms used for convergence tests in combination with
+    [`initial_condition_convergence_test`](@ref).
+    """
+    @inline function source_terms_convergence_test(u, x, t,
+                                                   equations::TrafficFlowLWREquations1D)
+        # Same settings as in `initial_condition`
+        c = 2.0
+        A = 1.0
+        L = 1
+        f = 1 / L
+        omega = 2 * pi * f
+        du = omega * cos(omega * (x[1] - t)) *
+             (-1 - 2 * equations.v_max * sin(omega * (x[1] - t)) - 3 * equations.v_max)
+
+        return SVector(du)
+    end
+
+    # Calculate 1D flux in for a single point
+    @inline function flux(u, orientation::Integer, equations::TrafficFlowLWREquations1D)
+        return SVector(equations.v_max * u[1] * (1.0 - u[1]))
+    end
+
+    # Calculate maximum wave speed for local Lax-Friedrichs-type dissipation
+    @inline function max_abs_speed_naive(u_ll, u_rr, orientation::Integer,
+                                         equations::TrafficFlowLWREquations1D)
+        λ_max = max(abs(equations.v_max * (1.0 - 2 * u_ll[1])),
+                    abs(equations.v_max * (1.0 - 2 * u_rr[1])))
+    end
+
+    # Calculate minimum and maximum wave speeds for HLL-type fluxes
+    @inline function min_max_speed_naive(u_ll, u_rr, orientation::Integer,
+                                         equations::TrafficFlowLWREquations1D)
+        jac_L = equations.v_max * (1.0 - 2 * u_ll[1])
+        jac_R = equations.v_max * (1.0 - 2 * u_rr[1])
+
+        λ_min = min(jac_L, jac_R)
+        λ_max = max(jac_L, jac_R)
+
+        return λ_min, λ_max
+    end
+
+    @inline function max_abs_speeds(u, equations::TrafficFlowLWREquations1D)
+        return (abs(equations.v_max * (1.0 - 2 * u[1])),)
+    end
+
+    # Convert conservative variables to primitive
+    @inline cons2prim(u, equations::TrafficFlowLWREquations1D) = u
+
+    # Convert conservative variables to entropy variables
+    @inline cons2entropy(u, equations::TrafficFlowLWREquations1D) = u
+
+    # Calculate entropy for a conservative state `cons`
+    @inline entropy(u::Real, ::TrafficFlowLWREquations1D) = 0.5 * u^2
+    @inline entropy(u, equations::TrafficFlowLWREquations1D) = entropy(u[1], equations)
+
+    # Calculate total energy for a conservative state `cons`
+    @inline energy_total(u::Real, ::TrafficFlowLWREquations1D) = 0.5 * u^2
+    @inline energy_total(u, equations::TrafficFlowLWREquations1D) = energy_total(u[1],
+                                                                                 equations)
 end # @muladd
diff --git a/test/test_tree_1d_traffic_flow_lwr.jl b/test/test_tree_1d_traffic_flow_lwr.jl
index 7a0aba8f938..165b3178157 100644
--- a/test/test_tree_1d_traffic_flow_lwr.jl
+++ b/test/test_tree_1d_traffic_flow_lwr.jl
@@ -11,9 +11,10 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_1d_dgsem")
 #! format: noindent
 
 @trixi_testset "elixir_traffic_flow_lwr_convergence.jl" begin
-    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_traffic_flow_lwr_convergence.jl"),
-                        l2 = [0.0008455067389588569], 
-                        linf = [0.004591951086623913])
+    @test_trixi_include(joinpath(EXAMPLES_DIR,
+                                 "elixir_traffic_flow_lwr_convergence.jl"),
+                        l2=[0.0008455067389588569],
+                        linf=[0.004591951086623913])
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)
     let
@@ -26,7 +27,7 @@ end
 
 @trixi_testset "elixir_traffic_flow_lwr_trafficjam.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_traffic_flow_lwr_trafficjam.jl"),
-                        l2 = [0.17591004141761846], linf = [0.5])
+                        l2=[0.17591004141761846], linf=[0.5])
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)
     let