Merge branch 'main' into NaNMath

.github/workflows/FormatCheck.yml

@@ -29,7 +29,7 @@ jobs:
         # TODO: Change the call below to
         #       format(".")
         run: |
-          julia  -e 'using Pkg; Pkg.add(PackageSpec(name = "JuliaFormatter"))'
+          julia  -e 'using Pkg; Pkg.add(PackageSpec(name = "JuliaFormatter", version="1.0.45"))'
           julia  -e 'using JuliaFormatter; format(["benchmark", "examples", "ext", "src", "test", "utils"])'
       - name: Format check
         run: |

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`
+- Added Lighthill-Whitham-Richards (LWR) traffic model
 
 ## Changes when updating to v0.6 from v0.5.x
 

Project.toml

@@ -1,7 +1,7 @@
 name = "Trixi"
 uuid = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb"
 authors = ["Michael Schlottke-Lakemper <[email protected]>", "Gregor Gassner <[email protected]>", "Hendrik Ranocha <[email protected]>", "Andrew R. Winters <[email protected]>", "Jesse Chan <[email protected]>"]
-version = "0.6.9-pre"
+version = "0.6.10-pre"
 
 [deps]
 CodeTracking = "da1fd8a2-8d9e-5ec2-8556-3022fb5608a2"
@@ -46,6 +46,7 @@ TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
 Triangulate = "f7e6ffb2-c36d-4f8f-a77e-16e897189344"
 TriplotBase = "981d1d27-644d-49a2-9326-4793e63143c3"
 TriplotRecipes = "808ab39a-a642-4abf-81ff-4cb34ebbffa3"
+TrixiBase = "9a0f1c46-06d5-4909-a5a3-ce25d3fa3284"
 UUIDs = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"
 
 [weakdeps]
@@ -98,6 +99,7 @@ TimerOutputs = "0.5.7"
 Triangulate = "2.0"
 TriplotBase = "0.1"
 TriplotRecipes = "0.1"
+TrixiBase = "0.1.1"
 UUIDs = "1.6"
 julia = "1.8"
 

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

benchmark/benchmarks.jl

@@ -2,6 +2,9 @@
 # readability
 #! format: off
 
+using Pkg
+Pkg.activate(@__DIR__)
+
 using BenchmarkTools
 using Trixi
 
@@ -47,13 +50,14 @@ end
 let
   SUITE["latency"] = BenchmarkGroup()
   SUITE["latency"]["default_example"] = @benchmarkable run(
-    `$(Base.julia_cmd()) -e 'using Trixi; trixi_include(default_example())'`) seconds=60
+    `$(Base.julia_cmd()) --project=$(@__DIR__) -e 'using Trixi; trixi_include(default_example())'`) seconds=60
   for polydeg in [3, 7]
     command = "using Trixi; trixi_include(joinpath(examples_dir(), \"tree_2d_dgsem\", \"elixir_advection_extended.jl\"), tspan=(0.0, 1.0e-10), polydeg=$(polydeg), save_restart=TrivialCallback(), save_solution=TrivialCallback())"
-    SUITE["latency"]["polydeg_$polydeg"] = @benchmarkable run($`$(Base.julia_cmd()) -e $command`) seconds=60
+    SUITE["latency"]["polydeg_$polydeg"] = @benchmarkable run(
+      $`$(Base.julia_cmd()) --project=$(@__DIR__) -e $command`) seconds=60
   end
   SUITE["latency"]["euler_2d"] = @benchmarkable run(
-    `$(Base.julia_cmd()) -e 'using Trixi; trixi_include(joinpath(examples_dir(), "tree_2d_dgsem", "elixir_euler_kelvin_helmholtz_instability.jl"), tspan=(0.0, 1.0e-10), save_restart=TrivialCallback(), save_solution=TrivialCallback())'`) seconds=60
+    `$(Base.julia_cmd()) --project=$(@__DIR__) -e 'using Trixi; trixi_include(joinpath(examples_dir(), "tree_2d_dgsem", "elixir_euler_kelvin_helmholtz_instability.jl"), tspan=(0.0, 1.0e-10), save_solution=TrivialCallback())'`) seconds=60
   SUITE["latency"]["mhd_2d"] = @benchmarkable run(
-    `$(Base.julia_cmd()) -e 'using Trixi; trixi_include(joinpath(examples_dir(), "tree_2d_dgsem", "elixir_mhd_blast_wave.jl"), tspan=(0.0, 1.0e-10), save_solution=TrivialCallback())'`) seconds=60
+    `$(Base.julia_cmd()) --project=$(@__DIR__) -e 'using Trixi; trixi_include(joinpath(examples_dir(), "tree_2d_dgsem", "elixir_mhd_blast_wave.jl"), tspan=(0.0, 1.0e-10), save_solution=TrivialCallback())'`) seconds=60
 end

benchmark/run_benchmarks.jl

@@ -1,3 +1,7 @@
+using Pkg
+Pkg.activate(@__DIR__)
+Pkg.develop(PackageSpec(path=dirname(@__DIR__)))
+Pkg.instantiate()
 
 using PkgBenchmark
 using Trixi

docs/Project.toml

@@ -10,6 +10,7 @@ OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 Trixi2Vtk = "bc1476a1-1ca6-4cc3-950b-c312b255ff95"
+TrixiBase = "9a0f1c46-06d5-4909-a5a3-ce25d3fa3284"
 
 [compat]
 CairoMakie = "0.6, 0.7, 0.8, 0.9, 0.10, 0.11"
@@ -23,3 +24,4 @@ OrdinaryDiffEq = "6.49.1"
 Plots = "1.9"
 Test = "1"
 Trixi2Vtk = "0.3"
+TrixiBase = "0.1.1"

docs/make.jl

@@ -8,6 +8,7 @@ end
 
 using Trixi
 using Trixi2Vtk
+using TrixiBase
 
 # Get Trixi.jl root directory
 trixi_root_dir = dirname(@__DIR__)
@@ -82,7 +83,7 @@ tutorials = create_tutorials(files)
 # Make documentation
 makedocs(
     # Specify modules for which docstrings should be shown
-    modules = [Trixi, Trixi2Vtk],
+    modules = [Trixi, TrixiBase, Trixi2Vtk],
     # Set sitename to Trixi.jl
     sitename = "Trixi.jl",
     # Provide additional formatting options
@@ -128,6 +129,7 @@ makedocs(
         "Troubleshooting and FAQ" => "troubleshooting.md",
         "Reference" => [
                         "Trixi.jl" => "reference-trixi.md",
+                        "TrixiBase.jl" => "reference-trixibase.md",
                         "Trixi2Vtk.jl" => "reference-trixi2vtk.md"
                        ],
         "Authors" => "authors.md",

docs/src/conventions.md

@@ -47,10 +47,12 @@ Trixi.jl is distributed with several examples in the form of elixirs, small
 Julia scripts containing everything to set up and run a simulation. Working
 interactively from the Julia REPL with these scripts can be quite convenient
 while for exploratory research and development of Trixi.jl. For example, you
-can use the convenience function [`trixi_include`](@ref) to `include` an elixir
-with some modified arguments. To enable this, it is helpful to use a consistent
-naming scheme in elixirs, since [`trixi_include`](@ref) can only perform simple
-replacements. Some standard variables names are
+can use the convenience function
+[`trixi_include`](@ref)
+to `include` an elixir with some modified arguments. To enable this, it is
+helpful to use a consistent naming scheme in elixirs, since
+[`trixi_include`](@ref)
+can only perform simple replacements. Some standard variables names are
 
 - `polydeg` for the polynomial degree of a solver
 - `surface_flux` for the numerical flux at surfaces

docs/src/reference-trixibase.md

@@ -0,0 +1,9 @@
+# TrixiBase.jl API
+
+```@meta
+CurrentModule = TrixiBase
+```
+
+```@autodocs
+Modules = [TrixiBase]
+```

examples/structured_1d_dgsem/elixir_traffic_flow_lwr_greenlight.jl

@@ -0,0 +1,80 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+
+equations = TrafficFlowLWREquations1D()
+
+solver = DGSEM(polydeg = 3, surface_flux = FluxHLL(min_max_speed_davis))
+
+coordinates_min = (-1.0,) # minimum coordinate
+coordinates_max = (1.0,) # maximum coordinate
+cells_per_dimension = (64,)
+
+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(coordinates_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_inflow,
+                       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(),
+            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,82 @@
+
+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
+
+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 = 1.0)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback, alive_callback,
+                        stepsize_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);
+
+summary_callback() # print the timer summary

src/Trixi.jl

@@ -70,6 +70,8 @@ using Triangulate: Triangulate, TriangulateIO, triangulate
 export TriangulateIO # for type parameter in DGMultiMesh
 using TriplotBase: TriplotBase
 using TriplotRecipes: DGTriPseudocolor
+@reexport using TrixiBase: trixi_include
+using TrixiBase: TrixiBase
 @reexport using SimpleUnPack: @unpack
 using SimpleUnPack: @pack!
 using DataStructures: BinaryHeap, FasterForward, extract_all!
@@ -135,7 +137,7 @@ include("callbacks_step/callbacks_step.jl")
 include("callbacks_stage/callbacks_stage.jl")
 include("semidiscretization/semidiscretization_euler_gravity.jl")
 
-# `trixi_include` and special elixirs such as `convergence_test`
+# Special elixirs such as `convergence_test`
 include("auxiliary/special_elixirs.jl")
 
 # Plot recipes and conversion functions to visualize results with Plots.jl
@@ -161,7 +163,8 @@ export AcousticPerturbationEquations2D,
        ShallowWaterTwoLayerEquations1D, ShallowWaterTwoLayerEquations2D,
        ShallowWaterEquationsQuasi1D,
        LinearizedEulerEquations2D,
-       PolytropicEulerEquations2D
+       PolytropicEulerEquations2D,
+       TrafficFlowLWREquations1D
 
 export LaplaceDiffusion1D, LaplaceDiffusion2D, LaplaceDiffusion3D,
        CompressibleNavierStokesDiffusion1D, CompressibleNavierStokesDiffusion2D,