From 3fef537218ce49c8618cf14164ee8743ea99b6f7 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Wed, 7 Feb 2024 16:31:32 +0100 Subject: [PATCH 01/18] Add classic LWR traffic flow model to Trixi --- .../elixir_traffic_flow_lwr_greenlight.jl | 87 +++++++++++++ .../elixir_traffic_flow_lwr_convergence.jl | 57 ++++++++ .../elixir_traffic_flow_lwr_trafficjam.jl | 83 ++++++++++++ src/Trixi.jl | 3 +- src/equations/equations.jl | 5 + src/equations/traffic_flow_lwr_1d.jl | 122 ++++++++++++++++++ test/test_structured_1d.jl | 15 +++ test/test_tree_1d.jl | 3 + test/test_tree_1d_traffic_flow_lwr.jl | 41 ++++++ 9 files changed, 415 insertions(+), 1 deletion(-) create mode 100644 examples/structured_1d_dgsem/elixir_traffic_flow_lwr_greenlight.jl create mode 100644 examples/tree_1d_dgsem/elixir_traffic_flow_lwr_convergence.jl create mode 100644 examples/tree_1d_dgsem/elixir_traffic_flow_lwr_trafficjam.jl create mode 100644 src/equations/traffic_flow_lwr_1d.jl create mode 100644 test/test_tree_1d_traffic_flow_lwr.jl diff --git a/examples/structured_1d_dgsem/elixir_traffic_flow_lwr_greenlight.jl b/examples/structured_1d_dgsem/elixir_traffic_flow_lwr_greenlight.jl new file mode 100644 index 00000000000..ea3102ab1e2 --- /dev/null +++ b/examples/structured_1d_dgsem/elixir_traffic_flow_lwr_greenlight.jl @@ -0,0 +1,87 @@ + +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 diff --git a/examples/tree_1d_dgsem/elixir_traffic_flow_lwr_convergence.jl b/examples/tree_1d_dgsem/elixir_traffic_flow_lwr_convergence.jl new file mode 100644 index 00000000000..3153f9977ab --- /dev/null +++ b/examples/tree_1d_dgsem/elixir_traffic_flow_lwr_convergence.jl @@ -0,0 +1,57 @@ + +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 diff --git a/examples/tree_1d_dgsem/elixir_traffic_flow_lwr_trafficjam.jl b/examples/tree_1d_dgsem/elixir_traffic_flow_lwr_trafficjam.jl new file mode 100644 index 00000000000..b2529026930 --- /dev/null +++ b/examples/tree_1d_dgsem/elixir_traffic_flow_lwr_trafficjam.jl @@ -0,0 +1,83 @@ + +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 diff --git a/src/Trixi.jl b/src/Trixi.jl index 8d74fbc9736..e161d9309bc 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -155,7 +155,8 @@ export AcousticPerturbationEquations2D, ShallowWaterTwoLayerEquations1D, ShallowWaterTwoLayerEquations2D, ShallowWaterEquationsQuasi1D, LinearizedEulerEquations2D, - PolytropicEulerEquations2D + PolytropicEulerEquations2D, + TrafficFlowLWREquations1D export LaplaceDiffusion1D, LaplaceDiffusion2D, LaplaceDiffusion3D, CompressibleNavierStokesDiffusion1D, CompressibleNavierStokesDiffusion2D, diff --git a/src/equations/equations.jl b/src/equations/equations.jl index c041bf117ba..62f2f1d2c27 100644 --- a/src/equations/equations.jl +++ b/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 diff --git a/src/equations/traffic_flow_lwr_1d.jl b/src/equations/traffic_flow_lwr_1d.jl new file mode 100644 index 00000000000..59d82f68700 --- /dev/null +++ b/src/equations/traffic_flow_lwr_1d.jl @@ -0,0 +1,122 @@ +# 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 -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_structured_1d.jl b/test/test_structured_1d.jl index f0eecfa9acd..fea06554c57 100644 --- a/test/test_structured_1d.jl +++ b/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 diff --git a/test/test_tree_1d.jl b/test/test_tree_1d.jl index 4654f6313f7..8b470278ffd 100644 --- a/test/test_tree_1d.jl +++ b/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 diff --git a/test/test_tree_1d_traffic_flow_lwr.jl b/test/test_tree_1d_traffic_flow_lwr.jl new file mode 100644 index 00000000000..7a0aba8f938 --- /dev/null +++ b/test/test_tree_1d_traffic_flow_lwr.jl @@ -0,0 +1,41 @@ +module TestExamples1DTrafficFlowLWR + +using Test +using Trixi + +include("test_trixi.jl") + +EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_1d_dgsem") + +@testset "Traffic-flow LWR" begin +#! 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]) + # 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 + +@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]) + # 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 + +end # module From 1597340271f0150f40f935ce4a0c0a376dba1043 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Wed, 7 Feb 2024 16:36:52 +0100 Subject: [PATCH 02/18] fmt --- .../elixir_traffic_flow_lwr_greenlight.jl | 42 ++-- .../elixir_traffic_flow_lwr_convergence.jl | 19 +- .../elixir_traffic_flow_lwr_trafficjam.jl | 40 ++-- src/equations/equations.jl | 4 +- src/equations/traffic_flow_lwr_1d.jl | 226 +++++++++--------- test/test_tree_1d_traffic_flow_lwr.jl | 9 +- 6 files changed, 163 insertions(+), 177 deletions(-) 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 From 04be105c143ca521ad05d48d8fd94c6c48f2c4a8 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Wed, 7 Feb 2024 16:49:03 +0100 Subject: [PATCH 03/18] shorten --- src/equations/traffic_flow_lwr_1d.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/equations/traffic_flow_lwr_1d.jl b/src/equations/traffic_flow_lwr_1d.jl index f2d4836ce90..967cafa2f72 100644 --- a/src/equations/traffic_flow_lwr_1d.jl +++ b/src/equations/traffic_flow_lwr_1d.jl @@ -66,7 +66,7 @@ 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) + (-1 - equations.v_max * (2 * sin(omega * (x[1] - t)) + 3)) return SVector(du) end From 8d9ad2e5ca0a567e10e118936fffda9bd45f2d1c Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Wed, 7 Feb 2024 22:34:49 +0100 Subject: [PATCH 04/18] Update examples/structured_1d_dgsem/elixir_traffic_flow_lwr_greenlight.jl --- .../structured_1d_dgsem/elixir_traffic_flow_lwr_greenlight.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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 34cc7f0fba0..b2b03e8d26a 100644 --- a/examples/structured_1d_dgsem/elixir_traffic_flow_lwr_greenlight.jl +++ b/examples/structured_1d_dgsem/elixir_traffic_flow_lwr_greenlight.jl @@ -46,7 +46,7 @@ function boundary_condition_outflow(u_inner, orientation, normal_direction, x, t return flux end -boundary_conditions = (x_neg = boundary_condition_outflow, +boundary_conditions = (x_neg = boundary_condition_inflow, x_pos = boundary_condition_outflow) initial_condition = initial_condition_greenlight From e82012da884ef5d3769725d8549e060f4c21069a Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Wed, 7 Feb 2024 22:49:14 +0100 Subject: [PATCH 05/18] Update examples/structured_1d_dgsem/elixir_traffic_flow_lwr_greenlight.jl --- .../structured_1d_dgsem/elixir_traffic_flow_lwr_greenlight.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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 b2b03e8d26a..bc6dc3ebf17 100644 --- a/examples/structured_1d_dgsem/elixir_traffic_flow_lwr_greenlight.jl +++ b/examples/structured_1d_dgsem/elixir_traffic_flow_lwr_greenlight.jl @@ -32,7 +32,7 @@ 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(coordinates_min, t, equations) end boundary_condition_inflow = BoundaryConditionDirichlet(inflow) From 3eb1b39c707ff972652df0ae238f000b413ba5ad Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Thu, 8 Feb 2024 08:52:39 +0100 Subject: [PATCH 06/18] rm IC const, fmt --- src/equations/traffic_flow_lwr_1d.jl | 211 +++++++++++++-------------- 1 file changed, 102 insertions(+), 109 deletions(-) diff --git a/src/equations/traffic_flow_lwr_1d.jl b/src/equations/traffic_flow_lwr_1d.jl index 967cafa2f72..891e3249966 100644 --- a/src/equations/traffic_flow_lwr_1d.jl +++ b/src/equations/traffic_flow_lwr_1d.jl @@ -3,114 +3,107 @@ # 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 +#! format: noindent + +@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 - - 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 + +varnames(::typeof(cons2cons), ::TrafficFlowLWREquations1D) = ("car-density",) +varnames(::typeof(cons2prim), ::TrafficFlowLWREquations1D) = ("car-density",) + +""" + 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 From 48ab702cda8a3cbc6930220acde7572cc0300c7a Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Thu, 8 Feb 2024 09:32:51 +0100 Subject: [PATCH 07/18] add davis wave speed estimate for upcoming change --- .../elixir_traffic_flow_lwr_greenlight.jl | 2 +- src/equations/traffic_flow_lwr_1d.jl | 5 +++++ 2 files changed, 6 insertions(+), 1 deletion(-) 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 bc6dc3ebf17..cdab67bee8e 100644 --- a/examples/structured_1d_dgsem/elixir_traffic_flow_lwr_greenlight.jl +++ b/examples/structured_1d_dgsem/elixir_traffic_flow_lwr_greenlight.jl @@ -7,7 +7,7 @@ using Trixi equations = TrafficFlowLWREquations1D() basis = LobattoLegendreBasis(3) -surface_flux = flux_hll +surface_flux = FluxHLL(min_max_speed_davis) solver = DGSEM(basis, surface_flux) coordinates_min = (-1.0,) # minimum coordinate diff --git a/src/equations/traffic_flow_lwr_1d.jl b/src/equations/traffic_flow_lwr_1d.jl index 891e3249966..04185b86543 100644 --- a/src/equations/traffic_flow_lwr_1d.jl +++ b/src/equations/traffic_flow_lwr_1d.jl @@ -88,6 +88,11 @@ end return λ_min, λ_max end +@inline function min_max_speed_davis(u_ll, u_rr, orientation::Integer, + equations::TrafficFlowLWREquations1D) + min_max_speed_naive(u_ll, u_rr, orientation, equations) +end + @inline function max_abs_speeds(u, equations::TrafficFlowLWREquations1D) return (abs(equations.v_max * (1.0 - 2 * u[1])),) end From d76c8e2accd03c44e06f59b7f248cfc18ec9406a Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Thu, 8 Feb 2024 09:37:24 +0100 Subject: [PATCH 08/18] news and md --- NEWS.md | 1 + README.md | 2 +- 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/NEWS.md b/NEWS.md index 02a723fca45..0fed69bf8f3 100644 --- a/NEWS.md +++ b/NEWS.md @@ -12,6 +12,7 @@ for human readability. - Different boundary conditions for quad/hex meshes in Abaqus format, even if not generated by HOHQMesh, can now be digested by Trixi in 2D and 3D. - Subcell (positivity) limiting support for nonlinear variables in 2D for `TreeMesh` +- LWR Traffic model ## Changes when updating to v0.6 from v0.5.x diff --git a/README.md b/README.md index c531ab4d1a4..71370d3478e 100644 --- a/README.md +++ b/README.md @@ -53,7 +53,7 @@ installation and postprocessing procedures. Its features include: * Hyperbolic diffusion equations for elliptic problems * Lattice-Boltzmann equations (D2Q9 and D3Q27 schemes) * Shallow water equations - * Several scalar conservation laws (e.g., linear advection, Burgers' equation) + * Several scalar conservation laws (e.g., linear advection, Burgers' equation, LWR traffic flow) * Multi-physics simulations * [Self-gravitating gas dynamics](https://github.com/trixi-framework/paper-self-gravitating-gas-dynamics) * Shared-memory parallelization via multithreading From 99a715263162b6a47fcb8afdcd20a59aa827f393 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Thu, 8 Feb 2024 09:43:06 +0100 Subject: [PATCH 09/18] Use euler --- .../tree_1d_dgsem/elixir_traffic_flow_lwr_trafficjam.jl | 6 +++--- test/test_tree_1d_traffic_flow_lwr.jl | 2 +- 2 files changed, 4 insertions(+), 4 deletions(-) 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 d92d91dc1b6..a5b0688da84 100644 --- a/examples/tree_1d_dgsem/elixir_traffic_flow_lwr_trafficjam.jl +++ b/examples/tree_1d_dgsem/elixir_traffic_flow_lwr_trafficjam.jl @@ -65,7 +65,7 @@ analysis_callback = AnalysisCallback(semi, interval = analysis_interval) alive_callback = AliveCallback(analysis_interval = analysis_interval) -stepsize_callback = StepsizeCallback(cfl = 2.0) +stepsize_callback = StepsizeCallback(cfl = 1.0) callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, @@ -74,8 +74,8 @@ callbacks = CallbackSet(summary_callback, ############################################################################### # run the simulation -sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), +sol = solve(ode, Euler(), 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 +summary_callback() # print the timer summary \ No newline at end of file diff --git a/test/test_tree_1d_traffic_flow_lwr.jl b/test/test_tree_1d_traffic_flow_lwr.jl index 165b3178157..54412e314b3 100644 --- a/test/test_tree_1d_traffic_flow_lwr.jl +++ b/test/test_tree_1d_traffic_flow_lwr.jl @@ -27,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.1761758135539748], linf=[0.5]) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) let From c71dea67546dfbe0ca83305e26201c5f7ff8ef1e Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Thu, 8 Feb 2024 10:23:26 +0100 Subject: [PATCH 10/18] fmt --- examples/tree_1d_dgsem/elixir_traffic_flow_lwr_trafficjam.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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 a5b0688da84..90af3bff6f9 100644 --- a/examples/tree_1d_dgsem/elixir_traffic_flow_lwr_trafficjam.jl +++ b/examples/tree_1d_dgsem/elixir_traffic_flow_lwr_trafficjam.jl @@ -78,4 +78,4 @@ sol = solve(ode, Euler(), 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 \ No newline at end of file +summary_callback() # print the timer summary From 4958216128184960fcd1b7ff58fbdd499bac2381 Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Thu, 8 Feb 2024 10:46:59 +0100 Subject: [PATCH 11/18] Update examples/structured_1d_dgsem/elixir_traffic_flow_lwr_greenlight.jl --- .../structured_1d_dgsem/elixir_traffic_flow_lwr_greenlight.jl | 1 - 1 file changed, 1 deletion(-) 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 cdab67bee8e..8d98eb6ddcb 100644 --- a/examples/structured_1d_dgsem/elixir_traffic_flow_lwr_greenlight.jl +++ b/examples/structured_1d_dgsem/elixir_traffic_flow_lwr_greenlight.jl @@ -14,7 +14,6 @@ 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) From ac73d1b5b0ec8981ec1641e933aa08a8f9a7a9c4 Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Thu, 8 Feb 2024 10:47:27 +0100 Subject: [PATCH 12/18] Update examples/tree_1d_dgsem/elixir_traffic_flow_lwr_trafficjam.jl --- examples/tree_1d_dgsem/elixir_traffic_flow_lwr_trafficjam.jl | 1 - 1 file changed, 1 deletion(-) 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 90af3bff6f9..7675e8e6dff 100644 --- a/examples/tree_1d_dgsem/elixir_traffic_flow_lwr_trafficjam.jl +++ b/examples/tree_1d_dgsem/elixir_traffic_flow_lwr_trafficjam.jl @@ -12,7 +12,6 @@ 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, From a3d70bc3366334320c39b823dabc3b17276ef5b8 Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Fri, 9 Feb 2024 00:00:04 +0100 Subject: [PATCH 13/18] Update NEWS.md Co-authored-by: Andrew Winters --- NEWS.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/NEWS.md b/NEWS.md index 0fed69bf8f3..06f8a80fbf6 100644 --- a/NEWS.md +++ b/NEWS.md @@ -12,7 +12,7 @@ for human readability. - Different boundary conditions for quad/hex meshes in Abaqus format, even if not generated by HOHQMesh, can now be digested by Trixi in 2D and 3D. - Subcell (positivity) limiting support for nonlinear variables in 2D for `TreeMesh` -- LWR Traffic model +- Added Lighthill-Whitham-Richards (LWR) traffic model ## Changes when updating to v0.6 from v0.5.x From 57488352d06d9e903abba8cff38ffe5ea901d975 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Fri, 9 Feb 2024 08:51:22 +0100 Subject: [PATCH 14/18] Andrew's suggestions --- examples/tree_1d_dgsem/elixir_traffic_flow_lwr_trafficjam.jl | 2 ++ src/equations/traffic_flow_lwr_1d.jl | 1 + 2 files changed, 3 insertions(+) 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 7675e8e6dff..d3a17b513fc 100644 --- a/examples/tree_1d_dgsem/elixir_traffic_flow_lwr_trafficjam.jl +++ b/examples/tree_1d_dgsem/elixir_traffic_flow_lwr_trafficjam.jl @@ -73,6 +73,8 @@ callbacks = CallbackSet(summary_callback, ############################################################################### # run the simulation +# Note: Be careful when increasing the polynomial degree and switching from first order finite volume +# to some actual DG method - in that case, you should also exchange the ODE solver. sol = solve(ode, Euler(), dt = 42, # solve needs some value here but it will be overwritten by the stepsize_callback save_everystep = false, callback = callbacks); diff --git a/src/equations/traffic_flow_lwr_1d.jl b/src/equations/traffic_flow_lwr_1d.jl index 04185b86543..44c51d4408c 100644 --- a/src/equations/traffic_flow_lwr_1d.jl +++ b/src/equations/traffic_flow_lwr_1d.jl @@ -9,6 +9,7 @@ TrafficFlowLWREquations1D The classic Lighthill-Witham Richards (LWR) model for 1D traffic flow. +The car density is denoted by $u$ and the maximum possible speed (e.g. due to speed limits) is $v_{\text{max}}$. ```math \partial_t u + v_{\text{max}} \partial_1 [u (1 - u)] = 0 ``` From 71e88a5e88941c354f2ee41e3f3fd4a0b32f4edf Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Tue, 13 Feb 2024 14:04:01 +0100 Subject: [PATCH 15/18] Apply suggestions from code review Co-authored-by: Hendrik Ranocha --- .../structured_1d_dgsem/elixir_traffic_flow_lwr_greenlight.jl | 4 +--- examples/tree_1d_dgsem/elixir_traffic_flow_lwr_convergence.jl | 4 +++- 2 files changed, 4 insertions(+), 4 deletions(-) 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 8d98eb6ddcb..e5badf14451 100644 --- a/examples/structured_1d_dgsem/elixir_traffic_flow_lwr_greenlight.jl +++ b/examples/structured_1d_dgsem/elixir_traffic_flow_lwr_greenlight.jl @@ -6,9 +6,7 @@ using Trixi equations = TrafficFlowLWREquations1D() -basis = LobattoLegendreBasis(3) -surface_flux = FluxHLL(min_max_speed_davis) -solver = DGSEM(basis, surface_flux) +solver = DGSEM(polydeg = 3, surface_flux = FluxHLL(min_max_speed_davis)) coordinates_min = (-1.0,) # minimum coordinate coordinates_max = (1.0,) # maximum coordinate 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 f887a198023..eb015f6fa9f 100644 --- a/examples/tree_1d_dgsem/elixir_traffic_flow_lwr_convergence.jl +++ b/examples/tree_1d_dgsem/elixir_traffic_flow_lwr_convergence.jl @@ -47,7 +47,9 @@ callbacks = CallbackSet(summary_callback, ############################################################################### # run the simulation -sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), +# Note: Be careful when increasing the polynomial degree and switching from first order finite volume +# to some actual DG method - in that case, you should also exchange the ODE solver. +sol = solve(ode, Euler(), dt = 42, # solve needs some value here but it will be overwritten by the stepsize_callback save_everystep = false, callback = callbacks); From 25a1b501b6f647d3fc181d923c2683791064b5ac Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Wed, 14 Feb 2024 09:19:11 +0100 Subject: [PATCH 16/18] add domain of u --- src/equations/traffic_flow_lwr_1d.jl | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/equations/traffic_flow_lwr_1d.jl b/src/equations/traffic_flow_lwr_1d.jl index 44c51d4408c..a4d2613a5c8 100644 --- a/src/equations/traffic_flow_lwr_1d.jl +++ b/src/equations/traffic_flow_lwr_1d.jl @@ -9,11 +9,12 @@ TrafficFlowLWREquations1D The classic Lighthill-Witham Richards (LWR) model for 1D traffic flow. -The car density is denoted by $u$ and the maximum possible speed (e.g. due to speed limits) is $v_{\text{max}}$. +The car density is denoted by $u \in [0, 1]$ and +the maximum possible speed (e.g. due to speed limits) is $v_{\text{max}}$. ```math \partial_t u + v_{\text{max}} \partial_1 [u (1 - u)] = 0 ``` -See e.g. Section 11.1 of +For more details 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 From ab8216f5c0b650d8cf41d08c29835d9bbb0a05a6 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Wed, 14 Feb 2024 17:58:36 +0100 Subject: [PATCH 17/18] back to carpenter kennedy --- examples/tree_1d_dgsem/elixir_traffic_flow_lwr_convergence.jl | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) 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 eb015f6fa9f..59258018f8c 100644 --- a/examples/tree_1d_dgsem/elixir_traffic_flow_lwr_convergence.jl +++ b/examples/tree_1d_dgsem/elixir_traffic_flow_lwr_convergence.jl @@ -47,9 +47,7 @@ callbacks = CallbackSet(summary_callback, ############################################################################### # run the simulation -# Note: Be careful when increasing the polynomial degree and switching from first order finite volume -# to some actual DG method - in that case, you should also exchange the ODE solver. -sol = solve(ode, Euler(), +sol = solve(ode, CarpenterKennedy2N54(), dt = 42, # solve needs some value here but it will be overwritten by the stepsize_callback save_everystep = false, callback = callbacks); From b0f4b0ff71e4fa78c289237426fdf0b51afc00da Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Thu, 15 Feb 2024 11:47:06 +0100 Subject: [PATCH 18/18] Update NEWS.md Co-authored-by: Andrew Winters --- NEWS.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/NEWS.md b/NEWS.md index 06f8a80fbf6..feccd1f9852 100644 --- a/NEWS.md +++ b/NEWS.md @@ -12,7 +12,7 @@ for human readability. - Different boundary conditions for quad/hex meshes in Abaqus format, even if not generated by HOHQMesh, can now be digested by Trixi in 2D and 3D. - Subcell (positivity) limiting support for nonlinear variables in 2D for `TreeMesh` -- Added Lighthill-Whitham-Richards (LWR) traffic model +- Added Lighthill-Whitham-Richards (LWR) traffic model ## Changes when updating to v0.6 from v0.5.x