Classic LWR traffic flow

examples/structured_1d_dgsem/elixir_traffic_flow_lwr_greenlight.jl

@@ -0,0 +1,83 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+
+equations = TrafficFlowLWREquations1D()
+
+basis = LobattoLegendreBasis(3)
+surface_flux = flux_hll
+solver = DGSEM(basis, surface_flux)
+
+coordinates_min = (-1.0,) # minimum coordinate
+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)
+
+# 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
+
+    return SVector(scalar)
+end
+
+###############################################################################
+# Specify non-periodic boundary conditions
+
+# Assume that there are always cars waiting at the left 
+function inflow(x, t, equations::TrafficFlowLWREquations1D)
+    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)
+
+    return flux
+end
+
+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)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 0.5)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 100
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+stepsize_callback = StepsizeCallback(cfl = 1.2)
+
+callbacks = CallbackSet(summary_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);
+
+summary_callback() # print the timer summary

examples/tree_1d_dgsem/elixir_traffic_flow_lwr_convergence.jl

@@ -0,0 +1,54 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+
+equations = TrafficFlowLWREquations1D()
+
+# Use first order finite volume to prevent oscillations at the shock
+solver = DGSEM(polydeg = 3, surface_flux = flux_hll)
+
+coordinates_min = 0.0 # minimum coordinate
+coordinates_max = 2.0 # maximum coordinate
+
+# Create a uniformly refined mesh with periodic boundaries
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level = 4,
+                n_cells_max = 30_000)
+
+###############################################################################
+# Specify non-periodic boundary conditions
+
+initial_condition = initial_condition_convergence_test
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
+                                    source_terms = source_terms_convergence_test)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 2.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)
+
+stepsize_callback = StepsizeCallback(cfl = 1.6)
+
+callbacks = CallbackSet(summary_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);
+
+summary_callback() # print the timer summary

examples/tree_1d_dgsem/elixir_traffic_flow_lwr_trafficjam.jl

@@ -0,0 +1,81 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+
+equations = TrafficFlowLWREquations1D()
+
+# Use first order finite volume to prevent oscillations at the shock
+solver = DGSEM(polydeg = 0, surface_flux = flux_lax_friedrichs)
+
+coordinates_min = -1.0 # minimum coordinate
+coordinates_max = 1.0 # maximum coordinate
+
+# Create a uniformly refined mesh with periodic boundaries
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level = 9,
+                n_cells_max = 30_000,
+                periodicity = false)
+
+# Example taken from http://www.clawpack.org/riemann_book/html/Traffic_flow.html#Example:-Traffic-jam                
+# 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
+
+    return SVector(scalar)
+end
+
+###############################################################################
+# Specify non-periodic boundary conditions
+
+function outflow(x, t, equations::TrafficFlowLWREquations1D)
+    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)
+
+    return flux
+end
+
+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)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 0.5)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 100
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+stepsize_callback = StepsizeCallback(cfl = 2.0)
+
+callbacks = CallbackSet(summary_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);
+
+summary_callback() # print the timer summary

src/Trixi.jl

@@ -155,7 +155,8 @@ export AcousticPerturbationEquations2D,
        ShallowWaterTwoLayerEquations1D, ShallowWaterTwoLayerEquations2D,
        ShallowWaterEquationsQuasi1D,
        LinearizedEulerEquations2D,
-       PolytropicEulerEquations2D
+       PolytropicEulerEquations2D,
+       TrafficFlowLWREquations1D
 
 export LaplaceDiffusion1D, LaplaceDiffusion2D, LaplaceDiffusion3D,
        CompressibleNavierStokesDiffusion1D, CompressibleNavierStokesDiffusion2D,

src/equations/equations.jl

@@ -507,4 +507,9 @@ include("linearized_euler_2d.jl")
 
 abstract type AbstractEquationsParabolic{NDIMS, NVARS, GradientVariables} <:
               AbstractEquations{NDIMS, NVARS} end
+
+# Lighthill-Witham-Richards (LWR) traffic flow model
+abstract type AbstractTrafficFlowLWREquations{NDIMS, NVARS} <:
+              AbstractEquations{NDIMS, NVARS} end
+include("traffic_flow_lwr_1d.jl")
 end # @muladd

src/equations/traffic_flow_lwr_1d.jl

@@ -0,0 +1,116 @@
+# By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
+# Since these FMAs can increase the performance of many numerical algorithms,
+# 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 - equations.v_max * (2 * sin(omega * (x[1] - t)) + 3))
+
+        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

test/test_structured_1d.jl

@@ -138,6 +138,21 @@ end
         @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
     end
 end
+
+@trixi_testset "elixir_traffic_flow_lwr_greenlight.jl" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR,
+                                 "elixir_traffic_flow_lwr_greenlight.jl"),
+                        l2=[0.2005523261652845],
+                        linf=[0.5052827913468407])
+    # Ensure that we do not have excessive memory allocations
+    # (e.g., from type instabilities)
+    let
+        t = sol.t[end]
+        u_ode = sol.u[end]
+        du_ode = similar(u_ode)
+        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
+    end
+end
 end
 
 # Clean up afterwards: delete Trixi.jl output directory

test/test_tree_1d.jl

@@ -47,6 +47,9 @@ isdir(outdir) && rm(outdir, recursive = true)
 
     # FDSBP methods on the TreeMesh
     include("test_tree_1d_fdsbp.jl")
+
+    # Traffic flow LWR
+    include("test_tree_1d_traffic_flow_lwr.jl")
 end
 
 # Coverage test for all initial conditions