diff --git a/.github/workflows/Downgrade.yml b/.github/workflows/Downgrade.yml
new file mode 100644
index 00000000000..c84b1026d1b
--- /dev/null
+++ b/.github/workflows/Downgrade.yml
@@ -0,0 +1,86 @@
+name: Downgrade
+
+on:
+  pull_request:
+    paths-ignore:
+      - 'AUTHORS.md'
+      - 'CITATION.bib'
+      - 'CONTRIBUTING.md'
+      - 'LICENSE.md'
+      - 'NEWS.md'
+      - 'README.md'
+      - '.zenodo.json'
+      - '.github/workflows/benchmark.yml'
+      - '.github/workflows/CompatHelper.yml'
+      - '.github/workflows/TagBot.yml'
+      - 'benchmark/**'
+      - 'docs/**'
+      - 'utils/**'
+  workflow_dispatch:
+
+# Cancel redundant CI tests automatically
+concurrency:
+  group: ${{ github.workflow }}-${{ github.ref }}
+  cancel-in-progress: true
+
+jobs:
+  downgrade_test:
+    if: "!contains(github.event.head_commit.message, 'skip ci')"
+    # We could also include the Julia version as in
+    # name: ${{ matrix.trixi_test }} - ${{ matrix.os }} - Julia ${{ matrix.version }} - ${{ matrix.arch }} - ${{ github.event_name }}
+    # to be more specific. However, that requires us updating the required CI tests whenever we update Julia.
+    name: Downgrade ${{ matrix.trixi_test }} - ${{ matrix.os }} - ${{ matrix.arch }} - ${{ github.event_name }}
+    runs-on: ${{ matrix.os }}
+    strategy:
+      fail-fast: false
+      matrix:
+        version:
+          - '1.9'
+          # - '~1.9.0-0' # including development versions
+          # - 'nightly'
+        os:
+          - ubuntu-latest
+        arch:
+          - x64
+        trixi_test:
+          # - tree_part1
+          # - tree_part2
+          # - tree_part3
+          # - tree_part4
+          # - tree_part5
+          # - tree_part6
+          # - structured
+          # - p4est_part1
+          # - p4est_part2
+          # - t8code_part1
+          # - unstructured_dgmulti
+          # - parabolic
+          # - paper_self_gravitating_gas_dynamics
+          # - misc_part1
+          # - misc_part2
+          # - performance_specializations_part1
+          # - performance_specializations_part2
+          # - mpi
+          - threaded
+    steps:
+      - uses: actions/checkout@v4
+      - uses: julia-actions/setup-julia@v1
+        with:
+          version: ${{ matrix.version }}
+          arch: ${{ matrix.arch }}
+      - run: julia -e 'using InteractiveUtils; versioninfo(verbose=true)'
+      - uses: julia-actions/cache@v1
+      - uses: julia-actions/julia-downgrade-compat@v1
+        with:
+          skip: LinearAlgebra,Printf,SparseArrays,DiffEqBase
+          projects: ., test
+      - uses: julia-actions/julia-buildpkg@v1
+        env:
+          PYTHON: ""
+      - name: Run tests without coverage
+        uses: julia-actions/julia-runtest@v1
+        with:
+          coverage: false
+        env:
+          PYTHON: ""
+          TRIXI_TEST: ${{ matrix.trixi_test }}
diff --git a/NEWS.md b/NEWS.md
index 38c02fb6866..ecc91581e9a 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -25,6 +25,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
 
diff --git a/Project.toml b/Project.toml
index af3e6b4d078..9bed045637a 100644
--- a/Project.toml
+++ b/Project.toml
@@ -62,17 +62,17 @@ DiffEqCallbacks = "2.25"
 Downloads = "1.6"
 EllipsisNotation = "1.0"
 FillArrays = "0.13.2, 1"
-ForwardDiff = "0.10.18"
-HDF5 = "0.14, 0.15, 0.16, 0.17"
+ForwardDiff = "0.10.24"
+HDF5 = "0.16.10, 0.17"
 IfElse = "0.1"
 LinearAlgebra = "1"
 LinearMaps = "2.7, 3.0"
-LoopVectorization = "0.12.118"
+LoopVectorization = "0.12.151"
 MPI = "0.20"
 Makie = "0.19, 0.20"
 MuladdMacro = "0.2.2"
-Octavian = "0.3.5"
-OffsetArrays = "1.3"
+Octavian = "0.3.21"
+OffsetArrays = "1.12"
 P4est = "0.4.9"
 Polyester = "0.7.5"
 PrecompileTools = "1.1"
@@ -81,19 +81,19 @@ RecipesBase = "1.1"
 Reexport = "1.0"
 Requires = "1.1"
 SciMLBase = "1.90, 2"
-Setfield = "0.8, 1"
+Setfield = "1"
 SimpleUnPack = "1.1"
 SparseArrays = "1"
-StartUpDG = "0.17"
-Static = "0.3, 0.4, 0.5, 0.6, 0.7, 0.8"
+StartUpDG = "0.17.7"
+Static = "0.8.7"
 StaticArrayInterface = "1.4"
-StaticArrays = "1"
-StrideArrays = "0.1.18"
-StructArrays = "0.6"
+StaticArrays = "1.5"
+StrideArrays = "0.1.26"
+StructArrays = "0.6.11"
 SummationByPartsOperators = "0.5.41"
 T8code = "0.4.3, 0.5"
 TimerOutputs = "0.5.7"
-Triangulate = "2.0"
+Triangulate = "2.2"
 TriplotBase = "0.1"
 TriplotRecipes = "0.1"
 TrixiBase = "0.1.1"
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
diff --git a/benchmark/benchmarks.jl b/benchmark/benchmarks.jl
index 15d1d96c05f..0d6fabcd4a9 100644
--- a/benchmark/benchmarks.jl
+++ b/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
diff --git a/benchmark/run_benchmarks.jl b/benchmark/run_benchmarks.jl
index 3a92a9ba700..7b8c25752f8 100644
--- a/benchmark/run_benchmarks.jl
+++ b/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
diff --git a/docs/literate/src/files/behind_the_scenes_simulation_setup.jl b/docs/literate/src/files/behind_the_scenes_simulation_setup.jl
new file mode 100644
index 00000000000..c93660e9bc1
--- /dev/null
+++ b/docs/literate/src/files/behind_the_scenes_simulation_setup.jl
@@ -0,0 +1,253 @@
+#src # Behind the scenes of a simulation setup
+
+# This tutorial will guide you through a simple Trixi.jl setup ("elixir"), giving an overview of
+# what happens in the background during the initialization of a simulation. While the setup
+# described herein does not cover all details, it involves relatively stable parts of Trixi.jl that
+# are unlikely to undergo significant changes in the near future. The goal is to clarify some of
+# the more fundamental, *technical* concepts that are applicable to a variety of
+# (also more complex) configurations.
+
+# Trixi.jl follows the [method of lines](http://www.scholarpedia.org/article/Method_of_lines) concept for solving partial differential equations (PDEs).
+# Firstly, the PDEs are reduced to a (potentially huge) system of
+# ordinary differential equations (ODEs) by discretizing the spatial derivatives. Subsequently,
+# these generated ODEs may be solved with methods available in OrdinaryDiffEq.jl or those specifically
+# implemented in Trixi.jl. The following steps elucidate the process of transitioning from PDEs to
+# ODEs within the framework of Trixi.jl.
+
+# ## Basic setup
+
+# Import essential libraries and specify an equation.
+
+using Trixi, OrdinaryDiffEq
+equations = LinearScalarAdvectionEquation2D((-0.2, 0.7))
+
+# Generate a spatial discretization using a [`TreeMesh`](@ref) with a pre-coarsened set of cells.
+
+coordinates_min = (-2.0, -2.0)
+coordinates_max = (2.0, 2.0)
+
+coarsening_patches = ((type = "box", coordinates_min = [0.0, -2.0],
+                       coordinates_max = [2.0, 0.0]),)
+
+mesh = TreeMesh(coordinates_min, coordinates_max, initial_refinement_level = 2,
+                n_cells_max = 30_000,
+                coarsening_patches = coarsening_patches)
+
+# The created `TreeMesh` looks like the following:
+
+# ![TreeMesh_example](https://github.com/trixi-framework/Trixi.jl/assets/119304909/d5ef76ee-8246-4730-a692-b472c06063a3)
+
+# Instantiate a [`DGSEM`](@ref) solver with a user-specified polynomial degree. The solver
+# will define `polydeg + 1` [Gauss-Lobatto nodes](https://en.wikipedia.org/wiki/Gaussian_quadrature#Gauss%E2%80%93Lobatto_rules) and their associated weights within
+# the reference interval ``[-1, 1]`` in each spatial direction. These nodes will be subsequently
+# used to approximate solutions on each leaf cell of the `TreeMesh`.
+
+solver = DGSEM(polydeg = 3)
+
+# Gauss-Lobatto nodes with `polydeg = 3`:
+
+# ![Gauss-Lobatto_nodes_example](https://github.com/trixi-framework/Trixi.jl/assets/119304909/1d894611-801e-4f75-bff0-d77ca1c672e5)
+
+# ## Overview of the [`SemidiscretizationHyperbolic`](@ref) type
+
+# At this stage, all necessary components for configuring the spatial discretization are in place.
+# The remaining task is to combine these components into a single structure that will be used
+# throughout the entire simulation process. This is where [`SemidiscretizationHyperbolic`](@ref)
+# comes into play.
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_convergence_test,
+                                    solver)
+
+# The constructor for the `SemidiscretizationHyperbolic` object calls numerous sub-functions to
+# perform the necessary initialization steps. A brief description of the key sub-functions is
+# provided below.
+
+
+# - `init_elements(leaf_cell_ids, mesh, equations, dg.basis, RealT, uEltype)`
+
+#   The fundamental elements for approximating the solution are the leaf
+#   cells. The solution is constructed as a polynomial of the degree specified in the `DGSEM`
+#   solver in each spatial direction on each leaf cell. This polynomial approximation is evaluated 
+#   at the Gauss-Lobatto nodes mentioned earlier. The `init_elements` function extracts
+#   these leaf cells from the `TreeMesh`, assigns them the label "elements", records their
+#   coordinates, and maps the Gauss-Lobatto nodes from the 1D interval ``[-1, 1]`` onto each coordinate axis
+#   of every element.
+
+
+#   ![elements_example](https://github.com/trixi-framework/Trixi.jl/assets/119304909/9f486670-b579-4e42-8697-439540c8bbb4)
+
+#   The visualization of elements with nodes shown here includes spaces between elements, which do
+#   not exist in reality. This spacing is included only for illustrative purposes to underscore the
+#   separation of elements and the independent projection of nodes onto each element.
+
+
+# - `init_interfaces(leaf_cell_ids, mesh, elements)`
+
+#   At this point, the elements with nodes have been defined; however, they lack the necessary
+#   communication functionality. This is crucial because the local solution polynomials on the
+#   elements are not independent of each other. Furthermore, nodes on the boundary of adjacent
+#   elements share the same spatial location, which requires a method to combine this into a
+#   meaningful solution.
+#   Here [Riemann solvers](https://en.wikipedia.org/wiki/Riemann_solver#Approximate_solvers)
+#   come into play which can handle the principal ambiguity of a multi-valued solution at the
+#   same spatial location.
+
+#   As demonstrated earlier, the elements can have varying sizes. Let us initially consider
+#   neighbors with equal size. For these elements, the `init_interfaces` function generates
+#   interfaces that store information about adjacent elements, their relative positions, and
+#   allocate containers for sharing solution data between neighbors during the solution process. 
+
+#   In our visualization, these interfaces would conceptually resemble tubes connecting the
+#   corresponding elements.
+
+#   ![interfaces_example](https://github.com/trixi-framework/Trixi.jl/assets/119304909/bc3b6b02-afbc-4371-aaf7-c7bdc5a6c540)
+
+
+# - `init_mortars(leaf_cell_ids, mesh, elements, dg.mortar)`
+
+#   Returning to the consideration of different sizes among adjacent elements, within the
+#   `TreeMesh`, adjacent leaf cells can vary in side length by a maximum factor of two. This
+#   implies that a large element has one neighbor of
+#   equal size with a connection through an interface, or two neighbors at half the size,
+#   requiring a connection through so called "mortars". In 3D, a large element would have
+#   four small neighbor elements.
+
+#   Mortars store information about the connected elements, their relative positions, and allocate
+#   containers for storing the solutions along the boundaries between these elements.
+
+#   Due to the differing sizes of adjacent elements, it is not feasible to directly map boundary
+#   nodes of adjacent elements. Therefore, the concept of mortars employs a mass-conserving
+#   interpolation function to map boundary nodes from a larger element to a smaller one.
+
+#   In our visualization, mortars are represented as branched tubes.
+
+#   ![mortars_example](https://github.com/trixi-framework/Trixi.jl/assets/119304909/43a95a60-3a31-4b1f-8724-14049e7a0481)
+
+
+# - `init_boundaries(leaf_cell_ids, mesh, elements)`
+
+#   In order to apply boundary conditions, it is necessary to identify the locations of the
+#   boundaries. Therefore, we initialize a "boundaries" object, which records the elements that
+#   contain boundaries, specifies which side of an element is a boundary, stores the coordinates
+#   of boundary nodes, and allocates containers for managing solutions at these boundaries.
+
+#   In our visualization, boundaries and their corresponding nodes are highlighted with green,
+#   semi-transparent lines.
+
+#   ![boundaries_example](https://github.com/trixi-framework/Trixi.jl/assets/119304909/21996b20-4a22-4dfb-b16a-e2c22c2f29fe)
+
+# All the structures mentioned earlier are collected as a cache of type `NamedTuple`. Subsequently,
+# an object of type `SemidiscretizationHyperbolic` is initialized using this cache, initial and
+# boundary conditions, equations, mesh and solver.
+
+# In conclusion, the primary purpose of a `SemidiscretizationHyperbolic` is to collect equations,
+# the geometric representation of the domain, and approximation instructions, creating specialized
+# structures to interconnect these components in a manner that enables their utilization for
+# the numerical solution of partial differential equations (PDEs).
+
+# As evident from the earlier description of `SemidiscretizationHyperbolic`, it comprises numerous
+# functions called subsequently. Without delving into details, the structure of the primary calls
+# are illustrated as follows:
+
+# ![SemidiscretizationHyperbolic_structure](https://github.com/trixi-framework/Trixi.jl/assets/119304909/8bf59422-0537-4d7a-9f13-d9b2253c19d7)
+
+# ## Overview of the [`semidiscretize`](@ref) function
+
+# At this stage, we have defined the equations and configured the domain's discretization. The
+# final step before solving is to select a suitable time span and apply the corresponding initial
+# conditions, which are already stored in the initialized `SemidiscretizationHyperbolic` object.
+
+# The purpose of the [`semidiscretize`](@ref) function is to wrap the semidiscretization as an
+# `ODEProblem` within the specified time interval. During this procedure the approximate solution
+#  is created at the given initial time via the specified `initial_condition` function from the 
+#  `SemidiscretizationHyperbolic` object. This `ODEProblem` can be subsequently passed to the
+# `solve` function from the [OrdinaryDiffEq.jl](https://github.com/SciML/OrdinaryDiffEq.jl) package
+# or to [`Trixi.solve`](@ref).
+
+ode = semidiscretize(semi, (0.0, 1.0));
+
+# The `semidiscretize` function involves a deep tree of subsequent calls, with the primary ones
+# explained below.
+
+
+# - `allocate_coefficients(mesh, equations, solver, cache)`
+
+#   To apply initial conditions, a data structure ("container") needs to be generated to store the
+#   initial values of the target variables for each node within each element.
+
+#   Since only one-dimensional `Array`s are `resize!`able in Julia, we use `Vector`s as an internal
+#   storage for the target variables and `resize!` them whenever needed, e.g. to change the number
+#   of elements. Then, during the solving process the same memory is reused by `unsafe_wrap`ping
+#   multi-dimensional `Array`s around the internal storage.
+
+# - `wrap_array(u_ode, semi)`
+
+#   As previously noted, `u_ode` is constructed as a 1D vector to ensure compatibility with
+#   OrdinaryDiffEq.jl. However, for internal use within Trixi.jl, identifying which part of the
+#   vector relates to specific variables, elements, or nodes can be challenging.
+
+#   This is why the `u_ode` vector is wrapped by the `wrap_array` function using `unsafe_wrap`
+#   to form a multidimensional array `u`. In this array, the first dimension corresponds to
+#   variables, followed by N dimensions corresponding to nodes for each of N space dimensions.
+#   The last dimension corresponds to the elements. 
+#   Consequently, navigation within this multidimensional array becomes noticeably easier.
+
+#   "Wrapping" in this context involves the creation of a reference to the same storage location
+#   but with an alternative structural representation. This approach enables the use of both
+#   instances `u` and `u_ode` as needed, so that changes are simultaneously reflected in both.
+#   This is possible because, from a storage perspective, they share the same stored data, while
+#   access to this data is provided in different ways.
+
+
+# - `compute_coefficients!(u, initial_conditions, t, mesh::DG, equations, solver, cache)`
+
+#   Now the variable `u`, intended to store solutions, has been allocated and wrapped, it is time
+#   to apply the initial conditions. The `compute_coefficients!` function calculates the initial
+#   conditions for each variable at every node within each element and properly stores them in the
+#   `u` array.
+
+# At this stage, the `semidiscretize` function has all the necessary components to initialize and
+# return an `ODEProblem` object, which will be used by the `solve` function to compute the
+# solution.
+
+# In summary, the internal workings of `semidiscretize` with brief descriptions can be presented
+# as follows.
+
+# ![semidiscretize_structure](https://github.com/trixi-framework/Trixi.jl/assets/119304909/491eddc4-aadb-4e29-8c76-a7c821d0674e)
+
+# ## Functions `solve` and `rhs!`
+
+# Once the `ODEProblem` object is initialized, the `solve` function and one of the ODE solvers from
+# the OrdinaryDiffEq.jl package can be utilized to compute an approximated solution using the
+# instructions contained in the `ODEProblem` object.
+
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), dt = 0.01,
+            save_everystep = false);
+
+# Since the `solve` function and the ODE solver have no knowledge
+# of a particular spatial discretization, it is necessary to define a
+# "right-hand-side function", `rhs!`, within Trixi.jl.
+
+# Trixi.jl includes a set of `rhs!` functions designed to compute `du`, i.e.,
+# ``\frac{\partial u}{\partial t}`` according to the structure
+# of the setup. These `rhs!` functions calculate interface, mortars, and boundary fluxes, in
+# addition to surface and volume integrals, in order to construct the `du` vector. This `du` vector
+# is then used by the time integration method to obtain the solution at the subsequent time step.
+# The `rhs!` function is called by time integration methods in each iteration of the solve loop
+# within OrdinaryDiffEq.jl, with arguments `du`, `u`, `semidiscretization`, and the current time.
+
+# Trixi.jl uses a two-levels approach for `rhs!` functions. The first level is limited to a
+# single function for each `semidiscretization` type, and its role is to redirect data to the
+# target `rhs!` for specific solver and mesh types. This target `rhs!` function is responsible
+# for calculating `du`.
+
+# Path from the `solve` function call to the appropriate `rhs!` function call:
+
+# ![rhs_structure](https://github.com/trixi-framework/Trixi.jl/assets/119304909/dbea9a0e-25a4-4afa-855e-01f1ad619982)
+
+# Computed solution:
+
+using Plots
+plot(sol)
+pd = PlotData2D(sol)
+plot!(getmesh(pd))
diff --git a/docs/literate/src/files/behind_the_scenes_simulation_setup_plots/Project.toml b/docs/literate/src/files/behind_the_scenes_simulation_setup_plots/Project.toml
new file mode 100644
index 00000000000..43aec5b7f54
--- /dev/null
+++ b/docs/literate/src/files/behind_the_scenes_simulation_setup_plots/Project.toml
@@ -0,0 +1,2 @@
+[deps]
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
diff --git a/docs/literate/src/files/behind_the_scenes_simulation_setup_plots/README.md b/docs/literate/src/files/behind_the_scenes_simulation_setup_plots/README.md
new file mode 100644
index 00000000000..011b5c75860
--- /dev/null
+++ b/docs/literate/src/files/behind_the_scenes_simulation_setup_plots/README.md
@@ -0,0 +1,15 @@
+# Plots for the tutorial "Behind the scenes of a simulation setup"
+
+To create all the images for the tutorial, execute the following command from the directory of this `README.md`:
+```julia
+pkg> activate .
+julia> include.(readdir("src"; join=true))
+```
+To create all images from a different directory, substitute `"src"` with the path to the `src` 
+folder. The resulting images will be generated in your current directory as PNG files.
+
+To generate a specific image, run the following command while replacing `"path/to/src"` and `"file_name"` with the appropriate values:
+```julia
+pkg> activate .
+julia> include(joinpath("path/to/src", "file_name"))
+```
\ No newline at end of file
diff --git a/docs/literate/src/files/behind_the_scenes_simulation_setup_plots/src/SemidiscretizationHyperbolic_structure_figure.jl b/docs/literate/src/files/behind_the_scenes_simulation_setup_plots/src/SemidiscretizationHyperbolic_structure_figure.jl
new file mode 100644
index 00000000000..cae7b19d470
--- /dev/null
+++ b/docs/literate/src/files/behind_the_scenes_simulation_setup_plots/src/SemidiscretizationHyperbolic_structure_figure.jl
@@ -0,0 +1,64 @@
+using Plots
+plot(Shape([(-2.3,4.5), (2.35,4.5), (2.35,2.5), (-2.3,2.5)]), linecolor="black", fillcolor="white", label=false,linewidth=2, size=(800,600), showaxis=false, grid=false, xlim=(-2.4,2.8), ylim=(-25,5.5))
+annotate!(2.3, 3.5, ("SemidiscretizationHyperbolic(mesh, equations, initial_conditions, solver; source_terms,
+boundary_conditions, RealT, uEltype, initial_cache)          ", 10, :black, :right))
+annotate!(-2.3, 1.5, ("creates and returns SemidiscretizationHyperbolic object, initialized using a mesh, equations,
+initial_conditions, boundary_conditions, source_terms, solver and cache", 9, :black, :left))
+plot!([-1.2,-1.2],[0.6,-2],arrow=true,color=:black,linewidth=2,label="")
+plot!([-1.2,-1.4],[0.6,-2],arrow=true,color=:black,linewidth=2,label="")
+plot!([-1.2,-1.],[0.6,-2],arrow=true,color=:black,linewidth=2,label="")
+annotate!(-1, -0.7, ("specialized for mesh
+and solver types", 9, :black, :left))
+plot!([1.25,1.25],[0.6,-2],arrow=true,color=:black,linewidth=2,label="")
+plot!([1.25,1.05],[0.6,-2],arrow=true,color=:black,linewidth=2,label="")
+plot!([1.25,1.45],[0.6,-2],arrow=true,color=:black,linewidth=2,label="")
+annotate!(1.48, -0.7, ("specialized for mesh
+and boundary_conditions
+types", 9, :black, :left))
+
+plot!(Shape([(-2.3,-2), (-0.1,-2), (-0.1,-4), (-2.3,-4)]), linecolor="black", fillcolor="white", label=false,linewidth=2)
+annotate!(-1.2, -3, ("create_cache(mesh::TreeMesh, equations,
+                   solver::Dg, RealT, uEltype)", 10, :black, :center))
+plot!([-2.22,-2.22],[-4,-22],arrow=false,color=:black,linewidth=2,label="")
+
+plot!(Shape([(-0.05,-2), (2.6,-2), (2.6,-4), (-0.05,-4)]), linecolor="black", fillcolor="white", label=false,linewidth=2)
+annotate!(1.27, -3, ("digest_boundary_conditions(boundary_conditions,
+                                             mesh, solver, cache)", 10, :black, :center))
+annotate!(2.6, -5, ("if necessary, converts passed boundary_conditions
+ into a suitable form for processing by Trixi.jl", 9, :black, :right))
+
+plot!(Shape([(-2,-6), (-0.55,-6), (-0.55,-7.1), (-2,-7.1)]), linecolor="black", fillcolor="white", label=false,linewidth=2)
+annotate!(-1.95, -6.5, ("local_leaf_cells(mesh.tree)", 10, :black, :left))
+annotate!(-2, -7.5, ("returns cells for which an element needs to be created (i.e. all leaf cells)", 9, :black, :left))
+plot!([-2.22,-2],[-6.5,-6.5],arrow=true,color=:black,linewidth=2,label="")
+
+plot!(Shape([(-2,-9), (1.73,-9), (1.73,-10.1), (-2,-10.1)]), linecolor="black", fillcolor="white", label=false,linewidth=2)
+annotate!(-1.95, -9.5, ("init_elements(leaf_cell_ids, mesh, equations, dg.basis, RealT, uEltype)", 10, :black, :left))
+annotate!(-2, -10.5, ("creates and initializes elements, projects Gauss-Lobatto basis onto each of them", 9, :black, :left))
+plot!([-2.22,-2],[-9.5,-9.5],arrow=true,color=:black,linewidth=2,label="")
+
+plot!(Shape([(-2,-12), (0.4,-12), (0.4,-13.1), (-2,-13.1)]), linecolor="black", fillcolor="white", label=false,linewidth=2)
+annotate!(-1.95, -12.5, ("init_interfaces(leaf_cell_ids, mesh, elements)", 10, :black, :left))
+annotate!(-2, -13.5, ("creates and initializes interfaces between each pair of adjacent elements of the same size", 9, :black, :left))
+plot!([-2.22,-2],[-12.5,-12.5],arrow=true,color=:black,linewidth=2,label="")
+
+plot!(Shape([(-2,-15), (0.5,-15), (0.5,-16.1), (-2,-16.1)]), linecolor="black", fillcolor="white", label=false,linewidth=2)
+annotate!(-1.95, -15.5, ("init_boundaries(leaf_cell_ids, mesh, elements)", 10, :black, :left))
+annotate!(-2, -17, ("creates and initializes boundaries, remembers each boundary element, as well as the coordinates of
+each boundary node", 9, :black, :left))
+plot!([-2.22,-2],[-15.5,-15.5],arrow=true,color=:black,linewidth=2,label="")
+
+plot!(Shape([(-1.6,-18), (1.3,-18), (1.3,-19.1), (-1.6,-19.1)]), linecolor="black", fillcolor="white", label=false,linewidth=2)
+annotate!(-1.55, -18.5, ("init_mortars(leaf_cell_ids, mesh, elements, dg.mortar)", 10, :black, :left))
+annotate!(-1.6, -20, ("creates and initializes mortars (type of interfaces) between each triple of adjacent coarsened
+and corresponding small elements", 9, :black, :left))
+plot!([-2.22,-1.6],[-18.5,-18.5],arrow=true,color=:black,linewidth=2,label="")
+annotate!(-2.15, -19, ("2D and 3D", 8, :black, :left))
+
+plot!(Shape([(-2,-21), (1.5,-21), (1.5,-23.1), (-2,-23.1)]), linecolor="black", fillcolor="white", label=false,linewidth=2)
+annotate!(-1.95, -22, ("create_cache(mesh, equations, dg.volume_integral, dg, uEltype)
+for 2D and 3D create_cache(mesh, equations, dg.mortar, uEltype)", 10, :black, :left))
+annotate!(-2, -23.5, ("add specialized parts of the cache required to compute the volume integral, etc.", 9, :black, :left))
+plot!([-2.22,-2],[-22,-22],arrow=true,color=:black,linewidth=2,label="")
+
+savefig("./SemidiscretizationHyperbolic")
\ No newline at end of file
diff --git a/docs/literate/src/files/behind_the_scenes_simulation_setup_plots/src/generate_boundary_figure.jl b/docs/literate/src/files/behind_the_scenes_simulation_setup_plots/src/generate_boundary_figure.jl
new file mode 100644
index 00000000000..14475d21339
--- /dev/null
+++ b/docs/literate/src/files/behind_the_scenes_simulation_setup_plots/src/generate_boundary_figure.jl
@@ -0,0 +1,190 @@
+using Plots
+
+function min(coordinates::Vector{Tuple{Float64, Float64}}, i)
+  min=coordinates[1][i]
+  for j in coordinates
+    if min>j[i]
+      min=j[i]
+    end
+  end
+  return min
+end
+
+function max(coordinates::Vector{Tuple{Float64, Float64}}, i)
+  max=coordinates[1][i]
+  for j in coordinates
+    if max<j[i]
+      max=j[i]
+    end
+  end
+  return max
+end
+
+function nodes_project_four(coordinates::Vector{Tuple{Float64, Float64}})
+  nodes_1d=[-1.0, -0.4472135954999579, 0.4472135954999579, 1.0]
+  nodes_2d=Array{Float64,2}(undef, 16,2)
+  for i in range(1,length(nodes_1d))
+    for j in range(1,length(nodes_1d))
+      nodes_2d[(i-1)*4+j,:]=[nodes_1d[i], nodes_1d[j]]
+    end
+  end
+  k_x = (max(coordinates,1)-min(coordinates,1))/2
+  k_y = (max(coordinates,2)-min(coordinates,2))/2
+  b_x = min(coordinates,1)
+  b_y = min(coordinates,2)
+  for i in range(1,Int(length(nodes_2d)/2))
+    nodes_2d[i,1]=k_x*(nodes_2d[i,1]+1)+b_x
+    nodes_2d[i,2]=k_y*(nodes_2d[i,2]+1)+b_y
+  end
+  return nodes_2d
+end
+
+
+import Base.+
++(a::Tuple,b::Tuple) = a .+ b
+
+# level 1
+plot(Shape([(2,0),(2,-2),(0,-2),(0,0)].+[(0.5,-0.75),(0.5,-0.75),(0.5,-0.75),(0.5,-0.75)]), linecolor="black", fillcolor="white", label="elements",legend_position=:right, linewidth=2, size=(800,600), showaxis=false, grid=false, x_lims=(-3,5.5), ylims=(-3.5,3))
+
+nodes_2d=nodes_project_four([(2,0),(2,-2),(0,-2),(0,0)].+[(0.5,-0.75),(0.5,-0.75),(0.5,-0.75),(0.5,-0.75)])
+scatter!(nodes_2d[:,1], nodes_2d[:,2], color = "red", label = "Gauss-Lobatto nodes", markersize = 5)
+
+# mortars
+plot!([-0.25,0.125], [-1,-1], linecolor="orange", fillcolor="white", label="mortars", linewidth=3)
+plot!([-0.25,0.125], [-2.5,-2.5], linecolor="orange", fillcolor="white", label=false, linewidth=3)
+plot!([0.125,0.125], [-1,-2.5], linecolor="orange", fillcolor="white", label=false, linewidth=3)
+plot!([0.125,0.5], [-1.75,-1.75], linecolor="orange", fillcolor="white", label=false, linewidth=3)
+plot!([0.75,0.75], [0,-0.325], linecolor="orange", fillcolor="white", label=false, linewidth=3)
+plot!([2.25,2.25], [0,-0.325], linecolor="orange", fillcolor="white", label=false, linewidth=3)
+plot!([2.25,0.75], [-0.325,-0.325], linecolor="orange", fillcolor="white", label=false, linewidth=3)
+plot!([1.5,1.5], [-0.75,-0.325], linecolor="orange", fillcolor="white", label=false, linewidth=3)
+
+plot!([2.5,2.5], [-0.75,-2.75], linecolor="green", label="boundaries with nodes", linewidth=15, opacity=0.3)
+plot!([.5,2.5], [-2.75,-2.75], linecolor="green", label=false, linewidth=15, opacity=0.3)
+
+# level 2
+
+# upper right
+# 1
+plot!(Shape([(2,2),(2,1),(1,1),(1,2)].+[(0.75,0.5),(0.75,0.5),(0.75,0.5),(0.75,0.5)]), linecolor="black", fillcolor="white", label=false, linewidth=2)
+
+nodes_2d=nodes_project_four([(2,2),(2,1),(1,1),(1,2)].+[(0.75,0.5),(0.75,0.5),(0.75,0.5),(0.75,0.5)])
+scatter!(nodes_2d[:,1], nodes_2d[:,2], color = "red", label = false, markersize = 5)
+
+plot!([1.25,1.75], [2,2], linecolor="blue", fillcolor="white", label="interfaces", linewidth=3)
+plot!([2.25,2.25], [1.5,1], linecolor="blue", fillcolor="white", label=false, linewidth=3)
+
+plot!([2.75,2.75], [2.5,1.5], linecolor="green", label=false, linewidth=15, opacity=0.3)
+plot!([1.75,2.75], [2.5,2.5], linecolor="green", label=false, linewidth=15, opacity=0.3)
+
+# 2
+plot!(Shape([(0,2),(1,2),(1,1),(0,1)].+[(0.25,0.5),(0.25,0.5),(0.25,0.5),(0.25,0.5)]), linecolor="black", fillcolor="white", label=false, linewidth=2)
+
+nodes_2d=nodes_project_four([(0,2),(1,2),(1,1),(0,1)].+[(0.25,0.5),(0.25,0.5),(0.25,0.5),(0.25,0.5)])
+scatter!(nodes_2d[:,1], nodes_2d[:,2], color = "red", label = false, markersize = 5)
+
+plot!([-0.25,0.25], [2,2], linecolor="blue", fillcolor="white", label=false, linewidth=3)
+plot!([0.75,0.75], [1.5,1], linecolor="blue", fillcolor="white", label=false, linewidth=3)
+
+plot!([.25,1.25], [2.5,2.5], linecolor="green", label=false, linewidth=15, opacity=0.3)
+
+# 3
+plot!(Shape([(0,0),(0,1),(1,1),(1,0)].+[(0.25,0.),(0.25,0.),(0.25,0.),(0.25,0.)]), linecolor="black", fillcolor="white", label=false, linewidth=2)
+
+nodes_2d=nodes_project_four([(0,0),(0,1),(1,1),(1,0)].+[(0.25,0.),(0.25,0.),(0.25,0.),(0.25,0.)])
+scatter!(nodes_2d[:,1], nodes_2d[:,2], color = "red", label = false, markersize = 5)
+
+plot!([1.25,1.75], [.5,.5], linecolor="blue", fillcolor="white", label=false, linewidth=3)
+plot!([-0.25,0.25], [.5,.5], linecolor="blue", fillcolor="white", label=false, linewidth=3)
+
+# 4
+plot!(Shape([(2,0),(2,1),(1,1),(1,0)].+[(0.75,0.),(0.75,0.),(0.75,0.),(0.75,0.)]), linecolor="black", fillcolor="white", label=false, linewidth=2)
+
+nodes_2d=nodes_project_four([(2,0),(2,1),(1,1),(1,0)].+[(0.75,0.),(0.75,0.),(0.75,0.),(0.75,0.)])
+scatter!(nodes_2d[:,1], nodes_2d[:,2], color = "red", label = false, markersize = 5)
+
+plot!([2.75,2.75], [1,0], linecolor="green", label=false, linewidth=15, opacity=0.3)
+
+# upper left
+# 1
+plot!(Shape([(-0,2),(-0,1),(-1,1),(-1,2)].+[(-0.25,0.5),(-0.25,0.5),(-0.25,0.5),(-0.25,0.5)]), linecolor="black", fillcolor="white", label=false, linewidth=2)
+
+nodes_2d=nodes_project_four([(-0,2),(-0,1),(-1,1),(-1,2)].+[(-0.25,0.5),(-0.25,0.5),(-0.25,0.5),(-0.25,0.5)])
+scatter!(nodes_2d[:,1], nodes_2d[:,2], color = "red", label = false, markersize = 5)
+
+plot!([-1.25,-1.75], [2,2], linecolor="blue", fillcolor="white", label=false, linewidth=3)
+plot!([-0.75,-0.75], [1,1.5], linecolor="blue", fillcolor="white", label=false, linewidth=3)
+
+plot!([-1.25,-0.25], [2.5,2.5], linecolor="green", label=false, linewidth=15, opacity=0.3)
+
+# 2
+plot!(Shape([(-2,2),(-1,2),(-1,1),(-2,1)].+[(-0.75,0.5),(-0.75,0.5),(-0.75,0.5),(-0.75,0.5)]), linecolor="black", fillcolor="white", label=false, linewidth=2)
+
+nodes_2d=nodes_project_four([(-2,2),(-1,2),(-1,1),(-2,1)].+[(-0.75,0.5),(-0.75,0.5),(-0.75,0.5),(-0.75,0.5)])
+scatter!(nodes_2d[:,1], nodes_2d[:,2], color = "red", label = false, markersize = 5)
+
+plot!([-2.25,-2.25], [1,1.5], linecolor="blue", fillcolor="white", label=false, linewidth=3)
+
+plot!([-1.75,-2.75], [2.5,2.5], linecolor="green", label=false, linewidth=15, opacity=0.3)
+plot!([-2.75,-2.75], [1.5,2.5], linecolor="green", label=false, linewidth=15, opacity=0.3)
+
+# 3
+plot!(Shape([(-2,0),(-2,1),(-1,1),(-1,0)].+[(-0.75,0.),(-0.75,0.),(-0.75,0.),(-0.75,0.)]), linecolor="black", fillcolor="white", label=false, linewidth=2)
+
+nodes_2d=nodes_project_four([(-2,0),(-2,1),(-1,1),(-1,0)].+[(-0.75,0.),(-0.75,0.),(-0.75,0.),(-0.75,0.)])
+scatter!(nodes_2d[:,1], nodes_2d[:,2], color = "red", label = false, markersize = 5)
+
+plot!([-1.75,-1.25], [.5,.5], linecolor="blue", fillcolor="white", label=false, linewidth=3)
+plot!([-2.25,-2.25], [-0.5,0], linecolor="blue", fillcolor="white", label=false, linewidth=3)
+
+plot!([-2.75,-2.75], [0,1], linecolor="green", label=false, linewidth=15, opacity=0.3)
+
+# 4
+plot!(Shape([(-0,0),(-1,0),(-1,1),(-0,1)].+[(-0.25,0.),(-0.25,0.),(-0.25,0.),(-0.25,0.)]), linecolor="black", fillcolor="white", label=false, linewidth=2)
+
+nodes_2d=nodes_project_four([(-0,0),(-1,0),(-1,1),(-0,1)].+[(-0.25,0.),(-0.25,0.),(-0.25,0.),(-0.25,0.)])
+scatter!(nodes_2d[:,1], nodes_2d[:,2], color = "red", label = false, markersize = 5)
+
+plot!([-0.75,-0.75], [-0.5,0], linecolor="blue", fillcolor="white", label=false, linewidth=3)
+
+# lower left
+# 1
+plot!(Shape([(-0,-0),(-0,-1),(-1,-1),(-1,-0)].+[(-0.25,-0.5),(-0.25,-0.5),(-0.25,-0.5),(-0.25,-0.5)]), linecolor="black", fillcolor="white", label=false, linewidth=2)
+
+nodes_2d=nodes_project_four([(-0,-0),(-0,-1),(-1,-1),(-1,-0)].+[(-0.25,-0.5),(-0.25,-0.5),(-0.25,-0.5),(-0.25,-0.5)])
+scatter!(nodes_2d[:,1], nodes_2d[:,2], color = "red", label = false, markersize = 5)
+
+plot!([-1.75,-1.25], [-1,-1], linecolor="blue", fillcolor="white", label=false, linewidth=3)
+plot!([-0.75,-0.75], [-1.5,-2], linecolor="blue", fillcolor="white", label=false, linewidth=3)
+
+
+# 2
+plot!(Shape([(-2,-1),(-2,0),(-1,-0),(-1,-1)].+[(-0.75,-0.5),(-0.75,-0.5),(-0.75,-0.5),(-0.75,-0.5)]), linecolor="black", fillcolor="white", label=false, linewidth=2)
+
+nodes_2d=nodes_project_four([(-2,-1),(-2,0),(-1,-0),(-1,-1)].+[(-0.75,-0.5),(-0.75,-0.5),(-0.75,-0.5),(-0.75,-0.5)])
+scatter!(nodes_2d[:,1], nodes_2d[:,2], color = "red", label = false, markersize = 5)
+
+plot!([-2.25,-2.25], [-1.5,-2], linecolor="blue", fillcolor="white", label=false, linewidth=3)
+
+plot!([-2.75,-2.75], [-0.5,-1.5], linecolor="green", label=false, linewidth=15, opacity=0.3)
+
+# 3
+plot!(Shape([(-2,-2),(-2,-1),(-1,-1),(-1,-2)].+[(-0.75,-1),(-0.75,-1),(-0.75,-1),(-0.75,-1)]), linecolor="black", fillcolor="white", label=false, linewidth=2)
+
+nodes_2d=nodes_project_four([(-2.,-2.),(-2.,-1.),(-1.,-1.),(-1.,-2.)].+[(-0.75,-1),(-0.75,-1),(-0.75,-1),(-0.75,-1)])
+scatter!(nodes_2d[:,1], nodes_2d[:,2], color = "red", label = false, markersize = 5)
+
+plot!([-1.75,-1.25], [-2.5,-2.5], linecolor="blue", fillcolor="white", label=false, linewidth=3)
+
+plot!([-2.75,-2.75], [-2,-3], linecolor="green", label=false, linewidth=15, opacity=0.3)
+plot!([-2.75,-1.75], [-3,-3], linecolor="green", label=false, linewidth=15, opacity=0.3)
+
+# 4
+plot!(Shape([(-1,-2),(-1,-1),(-0,-1),(-0,-2)].+[(-0.25,-1),(-0.25,-1),(-0.25,-1),(-0.25,-1)]), linecolor="black", fillcolor="white", label=false, linewidth=2)
+
+nodes_2d=nodes_project_four([(-1.,-2.),(-1.,-1.),(-0.,-1.),(-0.,-2.)].+[(-0.25,-1),(-0.25,-1),(-0.25,-1),(-0.25,-1)])
+scatter!(nodes_2d[:,1], nodes_2d[:,2], color = "red", label = false, markersize = 5)
+
+plot!([-0.25,-1.25], [-3,-3], linecolor="green", label=false, linewidth=15, opacity=0.3)
+
+savefig("./boundary")
\ No newline at end of file
diff --git a/docs/literate/src/files/behind_the_scenes_simulation_setup_plots/src/generate_elements_figure.jl b/docs/literate/src/files/behind_the_scenes_simulation_setup_plots/src/generate_elements_figure.jl
new file mode 100644
index 00000000000..6a98467a4e3
--- /dev/null
+++ b/docs/literate/src/files/behind_the_scenes_simulation_setup_plots/src/generate_elements_figure.jl
@@ -0,0 +1,117 @@
+using Plots
+
+function min(coordinates::Vector{Tuple{Float64, Float64}}, i)
+  min=coordinates[1][i]
+  for j in coordinates
+    if min>j[i]
+      min=j[i]
+    end
+  end
+  return min
+end
+
+function max(coordinates::Vector{Tuple{Float64, Float64}}, i)
+  max=coordinates[1][i]
+  for j in coordinates
+    if max<j[i]
+      max=j[i]
+    end
+  end
+  return max
+end
+
+function nodes_project_four(coordinates::Vector{Tuple{Float64, Float64}})
+  nodes_1d=[-1.0, -0.4472135954999579, 0.4472135954999579, 1.0]
+  nodes_2d=Array{Float64,2}(undef, 16,2)
+  for i in range(1,length(nodes_1d))
+    for j in range(1,length(nodes_1d))
+      nodes_2d[(i-1)*4+j,:]=[nodes_1d[i], nodes_1d[j]]
+    end
+  end
+  k_x = (max(coordinates,1)-min(coordinates,1))/2
+  k_y = (max(coordinates,2)-min(coordinates,2))/2
+  b_x = min(coordinates,1)
+  b_y = min(coordinates,2)
+  for i in range(1,Int(length(nodes_2d)/2))
+    nodes_2d[i,1]=k_x*(nodes_2d[i,1]+1)+b_x
+    nodes_2d[i,2]=k_y*(nodes_2d[i,2]+1)+b_y
+  end
+  return nodes_2d
+end
+
+
+import Base.+
++(a::Tuple,b::Tuple) = a .+ b
+
+# level 1
+plot(Shape([(2,0),(2,-2),(0,-2),(0,0)].+[(0.5,-0.75),(0.5,-0.75),(0.5,-0.75),(0.5,-0.75)]), linecolor="black", fillcolor="white", label="elements",legend_position=:right, linewidth=2, size=(800,600), showaxis=false, grid=false, x_lims=(-3,5.5), ylims=(-3.5,3))
+
+nodes_2d=nodes_project_four([(2,0),(2,-2),(0,-2),(0,0)].+[(0.5,-0.75),(0.5,-0.75),(0.5,-0.75),(0.5,-0.75)])
+scatter!(nodes_2d[:,1], nodes_2d[:,2], color = "red", label = "Gauss-Lobatto nodes", markersize = 5)
+
+# level 2
+
+# upper right
+plot!(Shape([(2,2),(2,1),(1,1),(1,2)].+[(0.75,0.5),(0.75,0.5),(0.75,0.5),(0.75,0.5)]), linecolor="black", fillcolor="white", label=false, linewidth=2)
+
+nodes_2d=nodes_project_four([(2,2),(2,1),(1,1),(1,2)].+[(0.75,0.5),(0.75,0.5),(0.75,0.5),(0.75,0.5)])
+scatter!(nodes_2d[:,1], nodes_2d[:,2], color = "red", label = false, markersize = 5)
+
+plot!(Shape([(0,2),(1,2),(1,1),(0,1)].+[(0.25,0.5),(0.25,0.5),(0.25,0.5),(0.25,0.5)]), linecolor="black", fillcolor="white", label=false, linewidth=2)
+
+nodes_2d=nodes_project_four([(0,2),(1,2),(1,1),(0,1)].+[(0.25,0.5),(0.25,0.5),(0.25,0.5),(0.25,0.5)])
+scatter!(nodes_2d[:,1], nodes_2d[:,2], color = "red", label = false, markersize = 5)
+
+plot!(Shape([(0,0),(0,1),(1,1),(1,0)].+[(0.25,0.),(0.25,0.),(0.25,0.),(0.25,0.)]), linecolor="black", fillcolor="white", label=false, linewidth=2)
+
+nodes_2d=nodes_project_four([(0,0),(0,1),(1,1),(1,0)].+[(0.25,0.),(0.25,0.),(0.25,0.),(0.25,0.)])
+scatter!(nodes_2d[:,1], nodes_2d[:,2], color = "red", label = false, markersize = 5)
+
+plot!(Shape([(2,0),(2,1),(1,1),(1,0)].+[(0.75,0.),(0.75,0.),(0.75,0.),(0.75,0.)]), linecolor="black", fillcolor="white", label=false, linewidth=2)
+
+nodes_2d=nodes_project_four([(2,0),(2,1),(1,1),(1,0)].+[(0.75,0.),(0.75,0.),(0.75,0.),(0.75,0.)])
+scatter!(nodes_2d[:,1], nodes_2d[:,2], color = "red", label = false, markersize = 5)
+
+# upper left
+plot!(Shape([(-0,2),(-0,1),(-1,1),(-1,2)].+[(-0.25,0.5),(-0.25,0.5),(-0.25,0.5),(-0.25,0.5)]), linecolor="black", fillcolor="white", label=false, linewidth=2)
+
+nodes_2d=nodes_project_four([(-0,2),(-0,1),(-1,1),(-1,2)].+[(-0.25,0.5),(-0.25,0.5),(-0.25,0.5),(-0.25,0.5)])
+scatter!(nodes_2d[:,1], nodes_2d[:,2], color = "red", label = false, markersize = 5)
+
+plot!(Shape([(-2,2),(-1,2),(-1,1),(-2,1)].+[(-0.75,0.5),(-0.75,0.5),(-0.75,0.5),(-0.75,0.5)]), linecolor="black", fillcolor="white", label=false, linewidth=2)
+
+nodes_2d=nodes_project_four([(-2,2),(-1,2),(-1,1),(-2,1)].+[(-0.75,0.5),(-0.75,0.5),(-0.75,0.5),(-0.75,0.5)])
+scatter!(nodes_2d[:,1], nodes_2d[:,2], color = "red", label = false, markersize = 5)
+
+plot!(Shape([(-2,0),(-2,1),(-1,1),(-1,0)].+[(-0.75,0.),(-0.75,0.),(-0.75,0.),(-0.75,0.)]), linecolor="black", fillcolor="white", label=false, linewidth=2)
+
+nodes_2d=nodes_project_four([(-2,0),(-2,1),(-1,1),(-1,0)].+[(-0.75,0.),(-0.75,0.),(-0.75,0.),(-0.75,0.)])
+scatter!(nodes_2d[:,1], nodes_2d[:,2], color = "red", label = false, markersize = 5)
+
+plot!(Shape([(-0,0),(-1,0),(-1,1),(-0,1)].+[(-0.25,0.),(-0.25,0.),(-0.25,0.),(-0.25,0.)]), linecolor="black", fillcolor="white", label=false, linewidth=2)
+
+nodes_2d=nodes_project_four([(-0,0),(-1,0),(-1,1),(-0,1)].+[(-0.25,0.),(-0.25,0.),(-0.25,0.),(-0.25,0.)])
+scatter!(nodes_2d[:,1], nodes_2d[:,2], color = "red", label = false, markersize = 5)
+
+# lower left
+plot!(Shape([(-0,-0),(-0,-1),(-1,-1),(-1,-0)].+[(-0.25,-0.5),(-0.25,-0.5),(-0.25,-0.5),(-0.25,-0.5)]), linecolor="black", fillcolor="white", label=false, linewidth=2)
+
+nodes_2d=nodes_project_four([(-0,-0),(-0,-1),(-1,-1),(-1,-0)].+[(-0.25,-0.5),(-0.25,-0.5),(-0.25,-0.5),(-0.25,-0.5)])
+scatter!(nodes_2d[:,1], nodes_2d[:,2], color = "red", label = false, markersize = 5)
+
+plot!(Shape([(-2,-1),(-2,0),(-1,-0),(-1,-1)].+[(-0.75,-0.5),(-0.75,-0.5),(-0.75,-0.5),(-0.75,-0.5)]), linecolor="black", fillcolor="white", label=false, linewidth=2)
+
+nodes_2d=nodes_project_four([(-2,-1),(-2,0),(-1,-0),(-1,-1)].+[(-0.75,-0.5),(-0.75,-0.5),(-0.75,-0.5),(-0.75,-0.5)])
+scatter!(nodes_2d[:,1], nodes_2d[:,2], color = "red", label = false, markersize = 5)
+
+plot!(Shape([(-2,-2),(-2,-1),(-1,-1),(-1,-2)].+[(-0.75,-1),(-0.75,-1),(-0.75,-1),(-0.75,-1)]), linecolor="black", fillcolor="white", label=false, linewidth=2)
+
+nodes_2d=nodes_project_four([(-2.,-2.),(-2.,-1.),(-1.,-1.),(-1.,-2.)].+[(-0.75,-1),(-0.75,-1),(-0.75,-1),(-0.75,-1)])
+scatter!(nodes_2d[:,1], nodes_2d[:,2], color = "red", label = false, markersize = 5)
+
+plot!(Shape([(-1,-2),(-1,-1),(-0,-1),(-0,-2)].+[(-0.25,-1),(-0.25,-1),(-0.25,-1),(-0.25,-1)]), linecolor="black", fillcolor="white", label=false, linewidth=2)
+
+nodes_2d=nodes_project_four([(-1.,-2.),(-1.,-1.),(-0.,-1.),(-0.,-2.)].+[(-0.25,-1),(-0.25,-1),(-0.25,-1),(-0.25,-1)])
+scatter!(nodes_2d[:,1], nodes_2d[:,2], color = "red", label = false, markersize = 5)
+
+savefig("./elements")
\ No newline at end of file
diff --git a/docs/literate/src/files/behind_the_scenes_simulation_setup_plots/src/generate_interfaces_figure.jl b/docs/literate/src/files/behind_the_scenes_simulation_setup_plots/src/generate_interfaces_figure.jl
new file mode 100644
index 00000000000..9d5e40722f8
--- /dev/null
+++ b/docs/literate/src/files/behind_the_scenes_simulation_setup_plots/src/generate_interfaces_figure.jl
@@ -0,0 +1,157 @@
+using Plots
+
+function min(coordinates::Vector{Tuple{Float64, Float64}}, i)
+  min=coordinates[1][i]
+  for j in coordinates
+    if min>j[i]
+      min=j[i]
+    end
+  end
+  return min
+end
+
+function max(coordinates::Vector{Tuple{Float64, Float64}}, i)
+  max=coordinates[1][i]
+  for j in coordinates
+    if max<j[i]
+      max=j[i]
+    end
+  end
+  return max
+end
+
+function nodes_project_four(coordinates::Vector{Tuple{Float64, Float64}})
+  nodes_1d=[-1.0, -0.4472135954999579, 0.4472135954999579, 1.0]
+  nodes_2d=Array{Float64,2}(undef, 16,2)
+  for i in range(1,length(nodes_1d))
+    for j in range(1,length(nodes_1d))
+      nodes_2d[(i-1)*4+j,:]=[nodes_1d[i], nodes_1d[j]]
+    end
+  end
+  k_x = (max(coordinates,1)-min(coordinates,1))/2
+  k_y = (max(coordinates,2)-min(coordinates,2))/2
+  b_x = min(coordinates,1)
+  b_y = min(coordinates,2)
+  for i in range(1,Int(length(nodes_2d)/2))
+    nodes_2d[i,1]=k_x*(nodes_2d[i,1]+1)+b_x
+    nodes_2d[i,2]=k_y*(nodes_2d[i,2]+1)+b_y
+  end
+  return nodes_2d
+end
+
+
+import Base.+
++(a::Tuple,b::Tuple) = a .+ b
+
+# level 1
+plot(Shape([(2,0),(2,-2),(0,-2),(0,0)].+[(0.5,-0.75),(0.5,-0.75),(0.5,-0.75),(0.5,-0.75)]), linecolor="black", fillcolor="white", label="elements",legend_position=:right, linewidth=2, size=(800,600), showaxis=false, grid=false, x_lims=(-3,5.5), ylims=(-3.5,3))
+
+nodes_2d=nodes_project_four([(2,0),(2,-2),(0,-2),(0,0)].+[(0.5,-0.75),(0.5,-0.75),(0.5,-0.75),(0.5,-0.75)])
+scatter!(nodes_2d[:,1], nodes_2d[:,2], color = "red", label = "Gauss-Lobatto nodes", markersize = 5)
+
+# level 2
+
+# upper right
+# 1
+plot!(Shape([(2,2),(2,1),(1,1),(1,2)].+[(0.75,0.5),(0.75,0.5),(0.75,0.5),(0.75,0.5)]), linecolor="black", fillcolor="white", label=false, linewidth=2)
+
+nodes_2d=nodes_project_four([(2,2),(2,1),(1,1),(1,2)].+[(0.75,0.5),(0.75,0.5),(0.75,0.5),(0.75,0.5)])
+scatter!(nodes_2d[:,1], nodes_2d[:,2], color = "red", label = false, markersize = 5)
+
+plot!([1.25,1.75], [2,2], linecolor="blue", fillcolor="white", label="interfaces", linewidth=3)
+plot!([2.25,2.25], [1.5,1], linecolor="blue", fillcolor="white", label=false, linewidth=3)
+
+# 2
+plot!(Shape([(0,2),(1,2),(1,1),(0,1)].+[(0.25,0.5),(0.25,0.5),(0.25,0.5),(0.25,0.5)]), linecolor="black", fillcolor="white", label=false, linewidth=2)
+
+nodes_2d=nodes_project_four([(0,2),(1,2),(1,1),(0,1)].+[(0.25,0.5),(0.25,0.5),(0.25,0.5),(0.25,0.5)])
+scatter!(nodes_2d[:,1], nodes_2d[:,2], color = "red", label = false, markersize = 5)
+
+plot!([-0.25,0.25], [2,2], linecolor="blue", fillcolor="white", label=false, linewidth=3)
+plot!([0.75,0.75], [1.5,1], linecolor="blue", fillcolor="white", label=false, linewidth=3)
+
+# 3
+plot!(Shape([(0,0),(0,1),(1,1),(1,0)].+[(0.25,0.),(0.25,0.),(0.25,0.),(0.25,0.)]), linecolor="black", fillcolor="white", label=false, linewidth=2)
+
+nodes_2d=nodes_project_four([(0,0),(0,1),(1,1),(1,0)].+[(0.25,0.),(0.25,0.),(0.25,0.),(0.25,0.)])
+scatter!(nodes_2d[:,1], nodes_2d[:,2], color = "red", label = false, markersize = 5)
+
+plot!([1.25,1.75], [.5,.5], linecolor="blue", fillcolor="white", label=false, linewidth=3)
+plot!([-0.25,0.25], [.5,.5], linecolor="blue", fillcolor="white", label=false, linewidth=3)
+
+# 4
+plot!(Shape([(2,0),(2,1),(1,1),(1,0)].+[(0.75,0.),(0.75,0.),(0.75,0.),(0.75,0.)]), linecolor="black", fillcolor="white", label=false, linewidth=2)
+
+nodes_2d=nodes_project_four([(2,0),(2,1),(1,1),(1,0)].+[(0.75,0.),(0.75,0.),(0.75,0.),(0.75,0.)])
+scatter!(nodes_2d[:,1], nodes_2d[:,2], color = "red", label = false, markersize = 5)
+
+# upper left
+# 1
+plot!(Shape([(-0,2),(-0,1),(-1,1),(-1,2)].+[(-0.25,0.5),(-0.25,0.5),(-0.25,0.5),(-0.25,0.5)]), linecolor="black", fillcolor="white", label=false, linewidth=2)
+
+nodes_2d=nodes_project_four([(-0,2),(-0,1),(-1,1),(-1,2)].+[(-0.25,0.5),(-0.25,0.5),(-0.25,0.5),(-0.25,0.5)])
+scatter!(nodes_2d[:,1], nodes_2d[:,2], color = "red", label = false, markersize = 5)
+
+plot!([-1.25,-1.75], [2,2], linecolor="blue", fillcolor="white", label=false, linewidth=3)
+plot!([-0.75,-0.75], [1,1.5], linecolor="blue", fillcolor="white", label=false, linewidth=3)
+
+# 2
+plot!(Shape([(-2,2),(-1,2),(-1,1),(-2,1)].+[(-0.75,0.5),(-0.75,0.5),(-0.75,0.5),(-0.75,0.5)]), linecolor="black", fillcolor="white", label=false, linewidth=2)
+
+nodes_2d=nodes_project_four([(-2,2),(-1,2),(-1,1),(-2,1)].+[(-0.75,0.5),(-0.75,0.5),(-0.75,0.5),(-0.75,0.5)])
+scatter!(nodes_2d[:,1], nodes_2d[:,2], color = "red", label = false, markersize = 5)
+
+plot!([-2.25,-2.25], [1,1.5], linecolor="blue", fillcolor="white", label=false, linewidth=3)
+
+
+# 3
+plot!(Shape([(-2,0),(-2,1),(-1,1),(-1,0)].+[(-0.75,0.),(-0.75,0.),(-0.75,0.),(-0.75,0.)]), linecolor="black", fillcolor="white", label=false, linewidth=2)
+
+nodes_2d=nodes_project_four([(-2,0),(-2,1),(-1,1),(-1,0)].+[(-0.75,0.),(-0.75,0.),(-0.75,0.),(-0.75,0.)])
+scatter!(nodes_2d[:,1], nodes_2d[:,2], color = "red", label = false, markersize = 5)
+
+plot!([-1.75,-1.25], [.5,.5], linecolor="blue", fillcolor="white", label=false, linewidth=3)
+plot!([-2.25,-2.25], [-0.5,0], linecolor="blue", fillcolor="white", label=false, linewidth=3)
+
+# 4
+plot!(Shape([(-0,0),(-1,0),(-1,1),(-0,1)].+[(-0.25,0.),(-0.25,0.),(-0.25,0.),(-0.25,0.)]), linecolor="black", fillcolor="white", label=false, linewidth=2)
+
+nodes_2d=nodes_project_four([(-0,0),(-1,0),(-1,1),(-0,1)].+[(-0.25,0.),(-0.25,0.),(-0.25,0.),(-0.25,0.)])
+scatter!(nodes_2d[:,1], nodes_2d[:,2], color = "red", label = false, markersize = 5)
+
+plot!([-0.75,-0.75], [-0.5,0], linecolor="blue", fillcolor="white", label=false, linewidth=3)
+
+# lower left
+# 1
+plot!(Shape([(-0,-0),(-0,-1),(-1,-1),(-1,-0)].+[(-0.25,-0.5),(-0.25,-0.5),(-0.25,-0.5),(-0.25,-0.5)]), linecolor="black", fillcolor="white", label=false, linewidth=2)
+
+nodes_2d=nodes_project_four([(-0,-0),(-0,-1),(-1,-1),(-1,-0)].+[(-0.25,-0.5),(-0.25,-0.5),(-0.25,-0.5),(-0.25,-0.5)])
+scatter!(nodes_2d[:,1], nodes_2d[:,2], color = "red", label = false, markersize = 5)
+
+plot!([-1.75,-1.25], [-1,-1], linecolor="blue", fillcolor="white", label=false, linewidth=3)
+plot!([-0.75,-0.75], [-1.5,-2], linecolor="blue", fillcolor="white", label=false, linewidth=3)
+
+
+# 2
+plot!(Shape([(-2,-1),(-2,0),(-1,-0),(-1,-1)].+[(-0.75,-0.5),(-0.75,-0.5),(-0.75,-0.5),(-0.75,-0.5)]), linecolor="black", fillcolor="white", label=false, linewidth=2)
+
+nodes_2d=nodes_project_four([(-2,-1),(-2,0),(-1,-0),(-1,-1)].+[(-0.75,-0.5),(-0.75,-0.5),(-0.75,-0.5),(-0.75,-0.5)])
+scatter!(nodes_2d[:,1], nodes_2d[:,2], color = "red", label = false, markersize = 5)
+
+plot!([-2.25,-2.25], [-1.5,-2], linecolor="blue", fillcolor="white", label=false, linewidth=3)
+
+# 3
+plot!(Shape([(-2,-2),(-2,-1),(-1,-1),(-1,-2)].+[(-0.75,-1),(-0.75,-1),(-0.75,-1),(-0.75,-1)]), linecolor="black", fillcolor="white", label=false, linewidth=2)
+
+nodes_2d=nodes_project_four([(-2.,-2.),(-2.,-1.),(-1.,-1.),(-1.,-2.)].+[(-0.75,-1),(-0.75,-1),(-0.75,-1),(-0.75,-1)])
+scatter!(nodes_2d[:,1], nodes_2d[:,2], color = "red", label = false, markersize = 5)
+
+plot!([-1.75,-1.25], [-2.5,-2.5], linecolor="blue", fillcolor="white", label=false, linewidth=3)
+
+# 4
+plot!(Shape([(-1,-2),(-1,-1),(-0,-1),(-0,-2)].+[(-0.25,-1),(-0.25,-1),(-0.25,-1),(-0.25,-1)]), linecolor="black", fillcolor="white", label=false, linewidth=2)
+
+nodes_2d=nodes_project_four([(-1.,-2.),(-1.,-1.),(-0.,-1.),(-0.,-2.)].+[(-0.25,-1),(-0.25,-1),(-0.25,-1),(-0.25,-1)])
+scatter!(nodes_2d[:,1], nodes_2d[:,2], color = "red", label = false, markersize = 5)
+
+savefig("./interfaces")
\ No newline at end of file
diff --git a/docs/literate/src/files/behind_the_scenes_simulation_setup_plots/src/generate_mortars_figure.jl b/docs/literate/src/files/behind_the_scenes_simulation_setup_plots/src/generate_mortars_figure.jl
new file mode 100644
index 00000000000..f48b410c42a
--- /dev/null
+++ b/docs/literate/src/files/behind_the_scenes_simulation_setup_plots/src/generate_mortars_figure.jl
@@ -0,0 +1,166 @@
+using Plots
+
+function min(coordinates::Vector{Tuple{Float64, Float64}}, i)
+  min=coordinates[1][i]
+  for j in coordinates
+    if min>j[i]
+      min=j[i]
+    end
+  end
+  return min
+end
+
+function max(coordinates::Vector{Tuple{Float64, Float64}}, i)
+  max=coordinates[1][i]
+  for j in coordinates
+    if max<j[i]
+      max=j[i]
+    end
+  end
+  return max
+end
+
+function nodes_project_four(coordinates::Vector{Tuple{Float64, Float64}})
+  nodes_1d=[-1.0, -0.4472135954999579, 0.4472135954999579, 1.0]
+  nodes_2d=Array{Float64,2}(undef, 16,2)
+  for i in range(1,length(nodes_1d))
+    for j in range(1,length(nodes_1d))
+      nodes_2d[(i-1)*4+j,:]=[nodes_1d[i], nodes_1d[j]]
+    end
+  end
+  k_x = (max(coordinates,1)-min(coordinates,1))/2
+  k_y = (max(coordinates,2)-min(coordinates,2))/2
+  b_x = min(coordinates,1)
+  b_y = min(coordinates,2)
+  for i in range(1,Int(length(nodes_2d)/2))
+    nodes_2d[i,1]=k_x*(nodes_2d[i,1]+1)+b_x
+    nodes_2d[i,2]=k_y*(nodes_2d[i,2]+1)+b_y
+  end
+  return nodes_2d
+end
+
+
+import Base.+
++(a::Tuple,b::Tuple) = a .+ b
+
+# level 1
+plot(Shape([(2,0),(2,-2),(0,-2),(0,0)].+[(0.5,-0.75),(0.5,-0.75),(0.5,-0.75),(0.5,-0.75)]), linecolor="black", fillcolor="white", label="elements",legend_position=:right, linewidth=2, size=(800,600), showaxis=false, grid=false, x_lims=(-3,5.5), ylims=(-3.5,3))
+
+nodes_2d=nodes_project_four([(2,0),(2,-2),(0,-2),(0,0)].+[(0.5,-0.75),(0.5,-0.75),(0.5,-0.75),(0.5,-0.75)])
+scatter!(nodes_2d[:,1], nodes_2d[:,2], color = "red", label = "Gauss-Lobatto nodes", markersize = 5)
+
+plot!([-0.25,0.125], [-1,-1], linecolor="orange", fillcolor="white", label="mortars", linewidth=3)
+plot!([-0.25,0.125], [-2.5,-2.5], linecolor="orange", fillcolor="white", label=false, linewidth=3)
+plot!([0.125,0.125], [-1,-2.5], linecolor="orange", fillcolor="white", label=false, linewidth=3)
+plot!([0.125,0.5], [-1.75,-1.75], linecolor="orange", fillcolor="white", label=false, linewidth=3)
+plot!([0.75,0.75], [0,-0.325], linecolor="orange", fillcolor="white", label=false, linewidth=3)
+plot!([2.25,2.25], [0,-0.325], linecolor="orange", fillcolor="white", label=false, linewidth=3)
+plot!([2.25,0.75], [-0.325,-0.325], linecolor="orange", fillcolor="white", label=false, linewidth=3)
+plot!([1.5,1.5], [-0.75,-0.325], linecolor="orange", fillcolor="white", label=false, linewidth=3)
+
+# level 2
+
+# upper right
+# 1
+plot!(Shape([(2,2),(2,1),(1,1),(1,2)].+[(0.75,0.5),(0.75,0.5),(0.75,0.5),(0.75,0.5)]), linecolor="black", fillcolor="white", label=false, linewidth=2)
+
+nodes_2d=nodes_project_four([(2,2),(2,1),(1,1),(1,2)].+[(0.75,0.5),(0.75,0.5),(0.75,0.5),(0.75,0.5)])
+scatter!(nodes_2d[:,1], nodes_2d[:,2], color = "red", label = false, markersize = 5)
+
+plot!([1.25,1.75], [2,2], linecolor="blue", fillcolor="white", label="interfaces", linewidth=3)
+plot!([2.25,2.25], [1.5,1], linecolor="blue", fillcolor="white", label=false, linewidth=3)
+
+# 2
+plot!(Shape([(0,2),(1,2),(1,1),(0,1)].+[(0.25,0.5),(0.25,0.5),(0.25,0.5),(0.25,0.5)]), linecolor="black", fillcolor="white", label=false, linewidth=2)
+
+nodes_2d=nodes_project_four([(0,2),(1,2),(1,1),(0,1)].+[(0.25,0.5),(0.25,0.5),(0.25,0.5),(0.25,0.5)])
+scatter!(nodes_2d[:,1], nodes_2d[:,2], color = "red", label = false, markersize = 5)
+
+plot!([-0.25,0.25], [2,2], linecolor="blue", fillcolor="white", label=false, linewidth=3)
+plot!([0.75,0.75], [1.5,1], linecolor="blue", fillcolor="white", label=false, linewidth=3)
+
+# 3
+plot!(Shape([(0,0),(0,1),(1,1),(1,0)].+[(0.25,0.),(0.25,0.),(0.25,0.),(0.25,0.)]), linecolor="black", fillcolor="white", label=false, linewidth=2)
+
+nodes_2d=nodes_project_four([(0,0),(0,1),(1,1),(1,0)].+[(0.25,0.),(0.25,0.),(0.25,0.),(0.25,0.)])
+scatter!(nodes_2d[:,1], nodes_2d[:,2], color = "red", label = false, markersize = 5)
+
+plot!([1.25,1.75], [.5,.5], linecolor="blue", fillcolor="white", label=false, linewidth=3)
+plot!([-0.25,0.25], [.5,.5], linecolor="blue", fillcolor="white", label=false, linewidth=3)
+
+# 4
+plot!(Shape([(2,0),(2,1),(1,1),(1,0)].+[(0.75,0.),(0.75,0.),(0.75,0.),(0.75,0.)]), linecolor="black", fillcolor="white", label=false, linewidth=2)
+
+nodes_2d=nodes_project_four([(2,0),(2,1),(1,1),(1,0)].+[(0.75,0.),(0.75,0.),(0.75,0.),(0.75,0.)])
+scatter!(nodes_2d[:,1], nodes_2d[:,2], color = "red", label = false, markersize = 5)
+
+# upper left
+# 1
+plot!(Shape([(-0,2),(-0,1),(-1,1),(-1,2)].+[(-0.25,0.5),(-0.25,0.5),(-0.25,0.5),(-0.25,0.5)]), linecolor="black", fillcolor="white", label=false, linewidth=2)
+
+nodes_2d=nodes_project_four([(-0,2),(-0,1),(-1,1),(-1,2)].+[(-0.25,0.5),(-0.25,0.5),(-0.25,0.5),(-0.25,0.5)])
+scatter!(nodes_2d[:,1], nodes_2d[:,2], color = "red", label = false, markersize = 5)
+
+plot!([-1.25,-1.75], [2,2], linecolor="blue", fillcolor="white", label=false, linewidth=3)
+plot!([-0.75,-0.75], [1,1.5], linecolor="blue", fillcolor="white", label=false, linewidth=3)
+
+# 2
+plot!(Shape([(-2,2),(-1,2),(-1,1),(-2,1)].+[(-0.75,0.5),(-0.75,0.5),(-0.75,0.5),(-0.75,0.5)]), linecolor="black", fillcolor="white", label=false, linewidth=2)
+
+nodes_2d=nodes_project_four([(-2,2),(-1,2),(-1,1),(-2,1)].+[(-0.75,0.5),(-0.75,0.5),(-0.75,0.5),(-0.75,0.5)])
+scatter!(nodes_2d[:,1], nodes_2d[:,2], color = "red", label = false, markersize = 5)
+
+plot!([-2.25,-2.25], [1,1.5], linecolor="blue", fillcolor="white", label=false, linewidth=3)
+
+
+# 3
+plot!(Shape([(-2,0),(-2,1),(-1,1),(-1,0)].+[(-0.75,0.),(-0.75,0.),(-0.75,0.),(-0.75,0.)]), linecolor="black", fillcolor="white", label=false, linewidth=2)
+
+nodes_2d=nodes_project_four([(-2,0),(-2,1),(-1,1),(-1,0)].+[(-0.75,0.),(-0.75,0.),(-0.75,0.),(-0.75,0.)])
+scatter!(nodes_2d[:,1], nodes_2d[:,2], color = "red", label = false, markersize = 5)
+
+plot!([-1.75,-1.25], [.5,.5], linecolor="blue", fillcolor="white", label=false, linewidth=3)
+plot!([-2.25,-2.25], [-0.5,0], linecolor="blue", fillcolor="white", label=false, linewidth=3)
+
+# 4
+plot!(Shape([(-0,0),(-1,0),(-1,1),(-0,1)].+[(-0.25,0.),(-0.25,0.),(-0.25,0.),(-0.25,0.)]), linecolor="black", fillcolor="white", label=false, linewidth=2)
+
+nodes_2d=nodes_project_four([(-0,0),(-1,0),(-1,1),(-0,1)].+[(-0.25,0.),(-0.25,0.),(-0.25,0.),(-0.25,0.)])
+scatter!(nodes_2d[:,1], nodes_2d[:,2], color = "red", label = false, markersize = 5)
+
+plot!([-0.75,-0.75], [-0.5,0], linecolor="blue", fillcolor="white", label=false, linewidth=3)
+
+# lower left
+# 1
+plot!(Shape([(-0,-0),(-0,-1),(-1,-1),(-1,-0)].+[(-0.25,-0.5),(-0.25,-0.5),(-0.25,-0.5),(-0.25,-0.5)]), linecolor="black", fillcolor="white", label=false, linewidth=2)
+
+nodes_2d=nodes_project_four([(-0,-0),(-0,-1),(-1,-1),(-1,-0)].+[(-0.25,-0.5),(-0.25,-0.5),(-0.25,-0.5),(-0.25,-0.5)])
+scatter!(nodes_2d[:,1], nodes_2d[:,2], color = "red", label = false, markersize = 5)
+
+plot!([-1.75,-1.25], [-1,-1], linecolor="blue", fillcolor="white", label=false, linewidth=3)
+plot!([-0.75,-0.75], [-1.5,-2], linecolor="blue", fillcolor="white", label=false, linewidth=3)
+
+
+# 2
+plot!(Shape([(-2,-1),(-2,0),(-1,-0),(-1,-1)].+[(-0.75,-0.5),(-0.75,-0.5),(-0.75,-0.5),(-0.75,-0.5)]), linecolor="black", fillcolor="white", label=false, linewidth=2)
+
+nodes_2d=nodes_project_four([(-2,-1),(-2,0),(-1,-0),(-1,-1)].+[(-0.75,-0.5),(-0.75,-0.5),(-0.75,-0.5),(-0.75,-0.5)])
+scatter!(nodes_2d[:,1], nodes_2d[:,2], color = "red", label = false, markersize = 5)
+
+plot!([-2.25,-2.25], [-1.5,-2], linecolor="blue", fillcolor="white", label=false, linewidth=3)
+
+# 3
+plot!(Shape([(-2,-2),(-2,-1),(-1,-1),(-1,-2)].+[(-0.75,-1),(-0.75,-1),(-0.75,-1),(-0.75,-1)]), linecolor="black", fillcolor="white", label=false, linewidth=2)
+
+nodes_2d=nodes_project_four([(-2.,-2.),(-2.,-1.),(-1.,-1.),(-1.,-2.)].+[(-0.75,-1),(-0.75,-1),(-0.75,-1),(-0.75,-1)])
+scatter!(nodes_2d[:,1], nodes_2d[:,2], color = "red", label = false, markersize = 5)
+
+plot!([-1.75,-1.25], [-2.5,-2.5], linecolor="blue", fillcolor="white", label=false, linewidth=3)
+
+# 4
+plot!(Shape([(-1,-2),(-1,-1),(-0,-1),(-0,-2)].+[(-0.25,-1),(-0.25,-1),(-0.25,-1),(-0.25,-1)]), linecolor="black", fillcolor="white", label=false, linewidth=2)
+
+nodes_2d=nodes_project_four([(-1.,-2.),(-1.,-1.),(-0.,-1.),(-0.,-2.)].+[(-0.25,-1),(-0.25,-1),(-0.25,-1),(-0.25,-1)])
+scatter!(nodes_2d[:,1], nodes_2d[:,2], color = "red", label = false, markersize = 5)
+
+savefig("./mortars")
\ No newline at end of file
diff --git a/docs/literate/src/files/behind_the_scenes_simulation_setup_plots/src/generate_nodes_figure.jl b/docs/literate/src/files/behind_the_scenes_simulation_setup_plots/src/generate_nodes_figure.jl
new file mode 100644
index 00000000000..b983563f47a
--- /dev/null
+++ b/docs/literate/src/files/behind_the_scenes_simulation_setup_plots/src/generate_nodes_figure.jl
@@ -0,0 +1,6 @@
+using Plots
+
+plot([-1,1],[0,0], linecolor="black", label="reference interval", size=(600, 200), ylim=(-0.5,0.5), yaxis=false, grid=false)
+scatter!([-1.0, -0.4472135954999579, 0.4472135954999579, 1.0],[0,0,0,0], color="red", label="Gauss-Lobatto nodes")
+
+savefig("./nodes")
\ No newline at end of file
diff --git a/docs/literate/src/files/behind_the_scenes_simulation_setup_plots/src/generate_treemesh_figure.jl b/docs/literate/src/files/behind_the_scenes_simulation_setup_plots/src/generate_treemesh_figure.jl
new file mode 100644
index 00000000000..8d2a8173a12
--- /dev/null
+++ b/docs/literate/src/files/behind_the_scenes_simulation_setup_plots/src/generate_treemesh_figure.jl
@@ -0,0 +1,26 @@
+using Plots
+# level 0
+plot(Shape([(2,2),(2,-2),(-2,-2),(-2,2)]), linecolor="black", fillcolor="white", label="0th level, root", linewidth=2, size=(600,600))
+# level 1
+plot!(Shape([(1.97,1.97),(1.97,0),(0,0),(0,1.97)]), linecolor="red", fillcolor="white", label="1st level", linewidth=2)
+plot!(Shape([(-1.97,1.97),(0,1.97),(0,0),(-1.97,0)]), linecolor="red", fillcolor="white", label=false, linewidth=2)
+plot!(Shape([(-1.97,0),(0,0),(0,-1.97),(-1.97,-1.97)]), linecolor="red", fillcolor="white", label=false, linewidth=2)
+plot!(Shape([(1.97,0),(1.97,-1.97),(0,-1.97),(0,0)]), linecolor="red", fillcolor="white", label=false, linewidth=2)
+# level 2
+# upper right
+plot!(Shape([(1.94,1.94),(1.94,1),(1,1),(1,1.94)]), linecolor="blue", fillcolor="white", label="2nd level", linewidth=2)
+plot!(Shape([(0.03,1.94),(1,1.94),(1,1),(0.03,1)]), linecolor="blue", fillcolor="white", label=false, linewidth=2)
+plot!(Shape([(0.03,0.03),(0.03,1),(1,1),(1,0.03)]), linecolor="blue", fillcolor="white", label=false, linewidth=2)
+plot!(Shape([(1.94,0.03),(1.94,1),(1,1),(1,0.03)]), linecolor="blue", fillcolor="white", label=false, linewidth=2)
+# upper left
+plot!(Shape([(-0.03,1.94),(-0.03,1),(-1,1),(-1,1.94)]), linecolor="blue", fillcolor="white", label=false, linewidth=2)
+plot!(Shape([(-1.94,1.94),(-1,1.94),(-1,1),(-1.94,1)]), linecolor="blue", fillcolor="white", label=false, linewidth=2)
+plot!(Shape([(-1.94,0.03),(-1.94,1),(-1,1),(-1,0.03)]), linecolor="blue", fillcolor="white", label=false, linewidth=2)
+plot!(Shape([(-0.03,0.03),(-1,0.03),(-1,1),(-0.03,1)]), linecolor="blue", fillcolor="white", label=false, linewidth=2)
+# lower left
+plot!(Shape([(-0.03,-0.03),(-0.03,-1),(-1,-1),(-1,-0.03)]), linecolor="blue", fillcolor="white", label=false, linewidth=2)
+plot!(Shape([(-1.94,-1),(-1.94,-0.03),(-1,-0.03),(-1,-1)]), linecolor="blue", fillcolor="white", label=false, linewidth=2)
+plot!(Shape([(-1.94,-1.94),(-1.94,-1),(-1,-1),(-1,-1.94)]), linecolor="blue", fillcolor="white", label=false, linewidth=2)
+plot!(Shape([(-1,-1.94),(-1,-1),(-0.03,-1),(-0.03,-1.94)]), linecolor="blue", fillcolor="white", label=false, linewidth=2)
+
+savefig("./treemesh")
\ No newline at end of file
diff --git a/docs/literate/src/files/behind_the_scenes_simulation_setup_plots/src/rhs_structure_figure.jl b/docs/literate/src/files/behind_the_scenes_simulation_setup_plots/src/rhs_structure_figure.jl
new file mode 100644
index 00000000000..19083dad68c
--- /dev/null
+++ b/docs/literate/src/files/behind_the_scenes_simulation_setup_plots/src/rhs_structure_figure.jl
@@ -0,0 +1,43 @@
+using Plots
+plot(Shape([(-0.4,0.5), (0.4,0.5), (0.4,-0.5), (-0.4,-0.5)]), linecolor="black", fillcolor="white", label=false,linewidth=2, size=(800,600), showaxis=false, grid=false, xlim=(-2,2), ylim=(-17.5,3))
+annotate!(0, 0, ("solve(...) call", 12, :black, :center))
+plot!([0,0],[-0.5,-2.5],arrow=true,color=:black,linewidth=2,label="")
+plot!(Shape([(-1.2,1), (1.2,1), (1.2,-1), (-1.2,-1)]), linecolor=:green, fillcolor=:transparent, label=false,linewidth=1, linestyle=:dash)
+annotate!(0.8, 0, ("Trixi.jl setup", 12, :green, :center))
+
+plot!(Shape([(-1.25,-2.5), (1.25,-2.5), (1.25,-3.5), (-1.25,-3.5)]), linecolor="black", fillcolor="white", label=false,linewidth=2)
+annotate!(0, -3, ("DiffEqBase.__solve(prob, alg, args...; kwargs...)", 12, :black, :center))
+plot!([0,0],[-3.5,-5.5],arrow=true,color=:black,linewidth=2,label="")
+plot!(Shape([(-1.9,-2), (1.9,-2), (1.9,-10), (-1.9,-10)]), linecolor=:red, fillcolor=:transparent, label=false,linewidth=1, linestyle=:dash)
+annotate!(1.4, -6, ("OrdinaryDiffEq.jl", 12, :red, :center))
+
+plot!(Shape([(-0.85,-5.5), (0.85,-5.5), (0.85,-6.5), (-0.85,-6.5)]), linecolor="black", fillcolor="white", label=false,linewidth=2)
+annotate!(0, -6, ("DiffEqBase.solve!(integrator)", 12, :black, :center))
+plot!([0,0],[-6.5,-8.5],arrow=true,color=:black,linewidth=2,label="")
+plot!([0,-0.5],[-6.5,-8.5],arrow=true,color=:black,linewidth=2,label="")
+plot!([0,0.5],[-6.5,-8.5],arrow=true,color=:black,linewidth=2,label="")
+annotate!(0.8, -7.5, ("Specialized for time \nintegration method", 9, :black, :center))
+
+plot!(Shape([(-1.2,-8.5), (1.2,-8.5), (1.2,-9.5), (-1.2,-9.5)]), linecolor="black", fillcolor="white", label=false,linewidth=2)
+annotate!(0, -9, ("perform_step!(integrator, integrator.cache)", 12, :black, :center))
+plot!([0,0],[-9.5,-11.5],arrow=true,color=:black,linewidth=2,label="")
+plot!([0,-0.5],[-9.5,-11.5],arrow=true,color=:black,linewidth=2,label="")
+plot!([0,0.5],[-9.5,-11.5],arrow=true,color=:black,linewidth=2,label="")
+annotate!(0.8, -10.5, ("Specialized for \nsemidiscretization", 9, :black, :center))
+
+plot!(Shape([(-0.95,-11.5), (0.95,-11.5), (0.95,-12.5), (-0.95,-12.5)]), linecolor="black", fillcolor="white", label=false,linewidth=2)
+annotate!(0, -12, ("Trixi.rhs!(du_ode, u_ode, semi, t)", 12, :black, :center))
+plot!([0,0],[-12.5,-14.5],arrow=true,color=:black,linewidth=2,label="")
+plot!([0,-0.5],[-12.5,-14.5],arrow=true,color=:black,linewidth=2,label="")
+plot!([0,0.5],[-12.5,-14.5],arrow=true,color=:black,linewidth=2,label="")
+annotate!(0.8, -13.5, ("Specialized for \nmesh and dg", 9, :black, :center))
+plot!(Shape([(-1.7,-11), (1.7,-11), (1.7,-17), (-1.7,-17)]), linecolor=:blue, fillcolor=:transparent, label=false,linewidth=1, linestyle=:dash)
+annotate!(1.5, -14, ("Trixi.jl", 12, :blue, :center))
+
+plot!(Shape([(-1.4,-14.5), (1.4,-14.5), (1.4,-16.5), (-1.4,-16.5)]), linecolor="black", fillcolor="white", label=false,linewidth=2)
+annotate!(0, -15.5, ("Trixi.rhs!(du, u, t, mesh, equations, initial_condition, \nboundary_conditions, source_terms, dg, cache)", 12, :black, :center))
+
+
+
+plot!([-2,2], [2,2], linecolor="white", fillcolor="white", label=false,linewidth=2)
+savefig("./rhs!")
\ No newline at end of file
diff --git a/docs/literate/src/files/behind_the_scenes_simulation_setup_plots/src/semidiscretize_structure_figure.jl b/docs/literate/src/files/behind_the_scenes_simulation_setup_plots/src/semidiscretize_structure_figure.jl
new file mode 100644
index 00000000000..52b0d573d7e
--- /dev/null
+++ b/docs/literate/src/files/behind_the_scenes_simulation_setup_plots/src/semidiscretize_structure_figure.jl
@@ -0,0 +1,51 @@
+using Plots
+plot(Shape([(-1.2,4), (1.2,4), (1.2,3), (-1.2,3)]), linecolor="black", fillcolor="white", label=false,linewidth=2, size=(800,600), showaxis=false, grid=false, xlim=(-2.6,2.6), ylim=(-21.5,4.5))
+annotate!(0, 3.5, ("semidiscretize(semi, tspan; reset_threads)", 10, :black, :center))
+annotate!(0, 2.5, ("creates and returns an ODEProblem object, initialized using rhs!, u0_ode, tspan and semi", 9, :black, :center))
+plot!([0,0],[2,0.5],arrow=true,color=:black,linewidth=2,label="")
+plot!([0,-0.2],[2,0.5],arrow=true,color=:black,linewidth=2,label="")
+plot!([0,0.2],[2,0.5],arrow=true,color=:black,linewidth=2,label="")
+annotate!(0.3,1.25, ("specialized for semi type", 9, :black, :left))
+
+plot!(Shape([(-1.6,0.5), (1.6,0.5), (1.6,-0.5), (-1.6,-0.5)]), linecolor="black", fillcolor="white", label=false,linewidth=2)
+annotate!(0, 0, ("compute_coefficients(t, semi::SemidiscretizationHyperbolic)", 10, :black, :center))
+plot!([0,0],[-0.5,-2],arrow=true,color=:black,linewidth=2,label="")
+
+plot!(Shape([(-1.3,-2), (1.3,-2), (1.3,-3), (-1.3,-3)]), linecolor="black", fillcolor="white", label=false,linewidth=2)
+annotate!(0, -2.5, ("compute_coefficients(initial_conditions, t, semi)", 10, :black, :center))
+plot!([0,0.2],[-3,-4.5],arrow=true,color=:black,linewidth=2,label="")
+plot!([0,-0.87],[-3,-9.5],arrow=true,color=:black,linewidth=2,label="")
+
+plot!(Shape([(0.1,-4.5), (2.5,-4.5), (2.5,-6.5), (0.1,-6.5)]), linecolor="black", fillcolor="white", label=false,linewidth=2)
+annotate!(0.15, -5.5, ("allocate_coefficients(mesh, equations, solver,
+                                 cache)", 10, :black, :left))
+annotate!(0.1, -8, ("initializes u_ode as a 1D zero-vector with a length
+depending on the number of variables, nodes,
+dimensions of mesh and elements, this vector is used
+by OrdinaryDiffEq.jl to keep a solution", 9, :black, :left))
+
+plot!(Shape([(-2.5,-9.5), (0.05,-9.5), (0.05,-11.5), (-2.5,-11.5)]), linecolor="black", fillcolor="white", label=false,linewidth=2)
+annotate!(-2.45, -10.5, ("compute_coefficients!(u_ode, initial_conditions,
+                                    t, semi)", 10, :black, :left))
+plot!([-2.4,-2.4],[-11.5,-18.5],arrow=false,color=:black,linewidth=2,label="")
+
+plot!(Shape([(-2.2,-13), (-0.8,-13), (-0.8,-14), (-2.2,-14)]), linecolor="black", fillcolor="white", label=false,linewidth=2)
+annotate!(-1.5, -13.5, ("wrap_array(u_ode, semi)", 10, :black, :center))
+annotate!(-2.2, -15.2, ("to simplify processing, it wraps the 1D array u_ode created in the allocate_coefficients function in
+a multidimensional one, where dimensions are responsible for variables, nodes on each axis and
+elements", 9, :black, :left))
+plot!([-2.4,-2.2],[-13.5,-13.5],arrow=true,color=:black,linewidth=2,label="")
+
+plot!(Shape([(-1.1,-17.5), (1.7,-17.5), (1.7,-19.5), (-1.1,-19.5)]), linecolor="black", fillcolor="white", label=false,linewidth=2)
+annotate!(-1.05, -18.5, ("compute_coefficients!(u, initial_conditions, t, mesh,
+                                    equations, solver, cache)", 10, :black, :left))
+annotate!(-1.1, -20.5, ("applies an initial conditions to each node of each element for each variable,
+saves in the wrapped u_ode", 9, :black, :left))
+plot!([-2.4,-1.1],[-18.5,-18.5],arrow=true,color=:black,linewidth=2,label="")
+plot!([-2.4,-1.1],[-18.5,-17.5],arrow=true,color=:black,linewidth=2,label="")
+plot!([-2.4,-1.1],[-18.5,-19.5],arrow=true,color=:black,linewidth=2,label="")
+annotate!(-2.4, -20.2, ("specialized for solver
+type and dimensionality
+of mesh", 9, :black, :left))
+
+savefig("./semidiscretize")
\ No newline at end of file
diff --git a/docs/literate/src/files/index.jl b/docs/literate/src/files/index.jl
index 1fc025d84da..6d5800c3b66 100644
--- a/docs/literate/src/files/index.jl
+++ b/docs/literate/src/files/index.jl
@@ -21,20 +21,29 @@
 # fundamental structure of a Trixi.jl setup, the visualization of results, and the development
 # process for Trixi.jl.
 
-# ### [2 Introduction to DG methods](@ref scalar_linear_advection_1d)
+# ### [2 Behind the scenes of a simulation setup](@ref behind_the_scenes_simulation_setup)
+#-
+# This tutorial will guide you through a simple Trixi.jl setup ("elixir"), giving an overview of
+# what happens in the background during the initialization of a simulation. While the setup
+# described herein does not cover all details, it involves relatively stable parts of Trixi.jl that
+# are unlikely to undergo significant changes in the near future. The goal is to clarify some of
+# the more fundamental, *technical* concepts that are applicable to a variety of
+# (also more complex) configurations.s
+
+# ### [3 Introduction to DG methods](@ref scalar_linear_advection_1d)
 #-
 # This tutorial gives an introduction to discontinuous Galerkin (DG) methods with the example of the
 # scalar linear advection equation in 1D. Starting with some theoretical explanations, we first implement
 # a raw version of a discontinuous Galerkin spectral element method (DGSEM). Then, we will show how
 # to use features of Trixi.jl to achieve the same result.
 
-# ### [3 DGSEM with flux differencing](@ref DGSEM_FluxDiff)
+# ### [4 DGSEM with flux differencing](@ref DGSEM_FluxDiff)
 #-
 # To improve stability often the flux differencing formulation of the DGSEM (split form) is used.
 # We want to present the idea and formulation on a basic 1D level. Then, we show how this formulation
 # can be implemented in Trixi.jl and analyse entropy conservation for two different flux combinations.
 
-# ### [4 Shock capturing with flux differencing and stage limiter](@ref shock_capturing)
+# ### [5 Shock capturing with flux differencing and stage limiter](@ref shock_capturing)
 #-
 # Using the flux differencing formulation, a simple procedure to capture shocks is a hybrid blending
 # of a high-order DG method and a low-order subcell finite volume (FV) method. We present the idea on a
@@ -42,20 +51,20 @@
 # explained and added to an exemplary simulation of the Sedov blast wave with the 2D compressible Euler
 # equations.
 
-# ### [5 Non-periodic boundary conditions](@ref non_periodic_boundaries)
+# ### [6 Non-periodic boundary conditions](@ref non_periodic_boundaries)
 #-
 # Thus far, all examples used periodic boundaries. In Trixi.jl, you can also set up a simulation with
 # non-periodic boundaries. This tutorial presents the implementation of the classical Dirichlet
 # boundary condition with a following example. Then, other non-periodic boundaries are mentioned.
 
-# ### [6 DG schemes via `DGMulti` solver](@ref DGMulti_1)
+# ### [7 DG schemes via `DGMulti` solver](@ref DGMulti_1)
 #-
 # This tutorial is about the more general DG solver [`DGMulti`](@ref), introduced [here](@ref DGMulti).
 # We are showing some examples for this solver, for instance with discretization nodes by Gauss or
 # triangular elements. Moreover, we present a simple way to include pre-defined triangulate meshes for
 # non-Cartesian domains using the package [StartUpDG.jl](https://github.com/jlchan/StartUpDG.jl).
 
-# ### [7 Other SBP schemes (FD, CGSEM) via `DGMulti` solver](@ref DGMulti_2)
+# ### [8 Other SBP schemes (FD, CGSEM) via `DGMulti` solver](@ref DGMulti_2)
 #-
 # Supplementary to the previous tutorial about DG schemes via the `DGMulti` solver we now present
 # the possibility for `DGMulti` to use other SBP schemes via the package
@@ -63,7 +72,7 @@
 # For instance, we show how to set up a finite differences (FD) scheme and a continuous Galerkin
 # (CGSEM) method.
 
-# ### [8 Upwind FD SBP schemes](@ref upwind_fdsbp)
+# ### [9 Upwind FD SBP schemes](@ref upwind_fdsbp)
 #-
 # General SBP schemes can not only be used via the [`DGMulti`](@ref) solver but
 # also with a general `DG` solver. In particular, upwind finite difference SBP
@@ -71,42 +80,42 @@
 # schemes in the `DGMulti` framework, the interface is based on the package
 # [SummationByPartsOperators.jl](https://github.com/ranocha/SummationByPartsOperators.jl).
 
-# ### [9 Adding a new scalar conservation law](@ref adding_new_scalar_equations)
+# ### [10 Adding a new scalar conservation law](@ref adding_new_scalar_equations)
 #-
 # This tutorial explains how to add a new physics model using the example of the cubic conservation
 # law. First, we define the equation using a `struct` `CubicEquation` and the physical flux. Then,
 # the corresponding standard setup in Trixi.jl (`mesh`, `solver`, `semi` and `ode`) is implemented
 # and the ODE problem is solved by OrdinaryDiffEq's `solve` method.
 
-# ### [10 Adding a non-conservative equation](@ref adding_nonconservative_equation)
+# ### [11 Adding a non-conservative equation](@ref adding_nonconservative_equation)
 #-
 # In this part, another physics model is implemented, the nonconservative linear advection equation.
 # We run two different simulations with different levels of refinement and compare the resulting errors.
 
-# ### [11 Parabolic terms](@ref parabolic_terms)
+# ### [12 Parabolic terms](@ref parabolic_terms)
 #-
 # This tutorial describes how parabolic terms are implemented in Trixi.jl, e.g.,
 # to solve the advection-diffusion equation.
 
-# ### [12 Adding new parabolic terms](@ref adding_new_parabolic_terms)
+# ### [13 Adding new parabolic terms](@ref adding_new_parabolic_terms)
 #-
 # This tutorial describes how new parabolic terms can be implemented using Trixi.jl.
 
-# ### [13 Adaptive mesh refinement](@ref adaptive_mesh_refinement)
+# ### [14 Adaptive mesh refinement](@ref adaptive_mesh_refinement)
 #-
 # Adaptive mesh refinement (AMR) helps to increase the accuracy in sensitive or turbolent regions while
 # not wasting resources for less interesting parts of the domain. This leads to much more efficient
 # simulations. This tutorial presents the implementation strategy of AMR in Trixi.jl, including the use of
 # different indicators and controllers.
 
-# ### [14 Structured mesh with curvilinear mapping](@ref structured_mesh_mapping)
+# ### [15 Structured mesh with curvilinear mapping](@ref structured_mesh_mapping)
 #-
 # In this tutorial, the use of Trixi.jl's structured curved mesh type [`StructuredMesh`](@ref) is explained.
 # We present the two basic option to initialize such a mesh. First, the curved domain boundaries
 # of a circular cylinder are set by explicit boundary functions. Then, a fully curved mesh is
 # defined by passing the transformation mapping.
 
-# ### [15 Unstructured meshes with HOHQMesh.jl](@ref hohqmesh_tutorial)
+# ### [16 Unstructured meshes with HOHQMesh.jl](@ref hohqmesh_tutorial)
 #-
 # The purpose of this tutorial is to demonstrate how to use the [`UnstructuredMesh2D`](@ref)
 # functionality of Trixi.jl. This begins by running and visualizing an available unstructured
@@ -115,26 +124,26 @@
 # software in the Trixi.jl ecosystem, and then run a simulation using Trixi.jl on said mesh.
 # In the end, the tutorial briefly explains how to simulate an example using AMR via `P4estMesh`.
 
-# ### [16 P4est mesh from gmsh](@ref p4est_from_gmsh)
+# ### [17 P4est mesh from gmsh](@ref p4est_from_gmsh)
 #-
 # This tutorial describes how to obtain a [`P4estMesh`](@ref) from an existing mesh generated
 # by [`gmsh`](https://gmsh.info/) or any other meshing software that can export to the Abaqus
 # input `.inp` format. The tutorial demonstrates how edges/faces can be associated with boundary conditions based on the physical nodesets.
 
-# ### [17 Explicit time stepping](@ref time_stepping)
+# ### [18 Explicit time stepping](@ref time_stepping)
 #-
 # This tutorial is about time integration using [OrdinaryDiffEq.jl](https://github.com/SciML/OrdinaryDiffEq.jl).
 # It explains how to use their algorithms and presents two types of time step choices - with error-based
 # and CFL-based adaptive step size control.
 
-# ### [18 Differentiable programming](@ref differentiable_programming)
+# ### [19 Differentiable programming](@ref differentiable_programming)
 #-
 # This part deals with some basic differentiable programming topics. For example, a Jacobian, its
 # eigenvalues and a curve of total energy (through the simulation) are calculated and plotted for
 # a few semidiscretizations. Moreover, we calculate an example for propagating errors with Measurement.jl
 # at the end.
 
-# ### [19 Custom semidiscretization](@ref custom_semidiscretization)
+# ### [20 Custom semidiscretization](@ref custom_semidiscretization)
 #-
 # This tutorial describes the [semidiscretiations](@ref overview-semidiscretizations) of Trixi.jl
 # and explains how to extend them for custom tasks.
diff --git a/docs/make.jl b/docs/make.jl
index 946b803b71e..8427c4049bf 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -55,6 +55,7 @@ files = [
         "Create first setup" => ("first_steps", "create_first_setup.jl"),
         "Changing Trixi.jl itself" => ("first_steps", "changing_trixi.jl"),
     ],
+    "Behind the scenes of a simulation setup" => "behind_the_scenes_simulation_setup.jl",
     # Topic: DG semidiscretizations
     "Introduction to DG methods" => "scalar_linear_advection_1d.jl",
     "DGSEM with flux differencing" => "DGSEM_FluxDiff.jl",
@@ -76,7 +77,7 @@ files = [
     # Topic: other stuff
     "Explicit time stepping" => "time_stepping.jl",
     "Differentiable programming" => "differentiable_programming.jl",
-    "Custom semidiscretizations" => "custom_semidiscretization.jl"
+    "Custom semidiscretizations" => "custom_semidiscretization.jl",
     ]
 tutorials = create_tutorials(files)
 
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..e5badf14451
--- /dev/null
+++ b/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
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..59258018f8c
--- /dev/null
+++ b/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
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..d3a17b513fc
--- /dev/null
+++ b/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
diff --git a/src/Trixi.jl b/src/Trixi.jl
index 866134a8fc4..efbd5eaa2f8 100644
--- a/src/Trixi.jl
+++ b/src/Trixi.jl
@@ -156,7 +156,8 @@ export AcousticPerturbationEquations2D,
        ShallowWaterEquations1D, ShallowWaterEquations2D,
        ShallowWaterEquationsQuasi1D,
        LinearizedEulerEquations2D,
-       PolytropicEulerEquations2D
+       PolytropicEulerEquations2D,
+       TrafficFlowLWREquations1D
 
 export LaplaceDiffusion1D, LaplaceDiffusion2D, LaplaceDiffusion3D,
        CompressibleNavierStokesDiffusion1D, CompressibleNavierStokesDiffusion2D,
diff --git a/src/equations/equations.jl b/src/equations/equations.jl
index d11ead8a2b8..8f476cf6f16 100644
--- a/src/equations/equations.jl
+++ b/src/equations/equations.jl
@@ -505,4 +505,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/polytropic_euler_2d.jl b/src/equations/polytropic_euler_2d.jl
index f5d2f7b0bad..e900fd64073 100644
--- a/src/equations/polytropic_euler_2d.jl
+++ b/src/equations/polytropic_euler_2d.jl
@@ -301,6 +301,46 @@ end
     return abs(v1) + c, abs(v2) + c
 end
 
+# Calculate maximum wave speed for local Lax-Friedrichs-type dissipation as the
+# maximum velocity magnitude plus the maximum speed of sound
+@inline function max_abs_speed_naive(u_ll, u_rr, orientation::Integer,
+                                     equations::PolytropicEulerEquations2D)
+    rho_ll, v1_ll, v2_ll = cons2prim(u_ll, equations)
+    rho_rr, v1_rr, v2_rr = cons2prim(u_rr, equations)
+
+    # Get the velocity value in the appropriate direction
+    if orientation == 1
+        v_ll = v1_ll
+        v_rr = v1_rr
+    else # orientation == 2
+        v_ll = v2_ll
+        v_rr = v2_rr
+    end
+    # Calculate sound speeds (we have p = kappa * rho^gamma)
+    c_ll = sqrt(equations.gamma * equations.kappa * rho_ll^(equations.gamma - 1))
+    c_rr = sqrt(equations.gamma * equations.kappa * rho_rr^(equations.gamma - 1))
+
+    λ_max = max(abs(v_ll), abs(v_rr)) + max(c_ll, c_rr)
+end
+
+@inline function max_abs_speed_naive(u_ll, u_rr, normal_direction::AbstractVector,
+                                     equations::PolytropicEulerEquations2D)
+    rho_ll, v1_ll, v2_ll = cons2prim(u_ll, equations)
+    rho_rr, v1_rr, v2_rr = cons2prim(u_rr, equations)
+
+    # Calculate normal velocities and sound speed (we have p = kappa * rho^gamma)
+    # left
+    v_ll = (v1_ll * normal_direction[1] +
+            v2_ll * normal_direction[2])
+    c_ll = sqrt(equations.gamma * equations.kappa * rho_ll^(equations.gamma - 1))
+    # right
+    v_rr = (v1_rr * normal_direction[1] +
+            v2_rr * normal_direction[2])
+    c_rr = sqrt(equations.gamma * equations.kappa * rho_rr^(equations.gamma - 1))
+
+    return max(abs(v_ll), abs(v_rr)) + max(c_ll, c_rr) * norm(normal_direction)
+end
+
 # Convert conservative variables to primitive
 @inline function cons2prim(u, equations::PolytropicEulerEquations2D)
     rho, rho_v1, rho_v2 = u
diff --git a/src/equations/traffic_flow_lwr_1d.jl b/src/equations/traffic_flow_lwr_1d.jl
new file mode 100644
index 00000000000..a4d2613a5c8
--- /dev/null
+++ b/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
+#! format: noindent
+
+@doc raw"""
+    TrafficFlowLWREquations1D
+
+The classic Lighthill-Witham Richards (LWR) model for 1D traffic flow.
+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
+```
+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
+"""
+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_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 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
+
+# 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/Project.toml b/test/Project.toml
index ecae0ac0900..a376c2805ea 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -13,13 +13,13 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [compat]
 Aqua = "0.8"
-CairoMakie = "0.6, 0.7, 0.8, 0.9, 0.10"
+CairoMakie = "0.10"
 Downloads = "1"
-ForwardDiff = "0.10"
+ForwardDiff = "0.10.24"
 LinearAlgebra = "1"
 MPI = "0.20"
 OrdinaryDiffEq = "6.49.1"
-Plots = "1.16"
+Plots = "1.19"
 Printf = "1"
 Random = "1"
 Test = "1"
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 66385fa6a2d..4a25a51a45e 100644
--- a/test/test_tree_1d.jl
+++ b/test/test_tree_1d.jl
@@ -45,6 +45,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..54412e314b3
--- /dev/null
+++ b/test/test_tree_1d_traffic_flow_lwr.jl
@@ -0,0 +1,42 @@
+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.1761758135539748], 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
diff --git a/test/test_unit.jl b/test/test_unit.jl
index 9614960f22e..1907a281718 100644
--- a/test/test_unit.jl
+++ b/test/test_unit.jl
@@ -880,6 +880,30 @@ end
     end
 end
 
+@timed_testset "Consistency check for Lax-Friedrich flux: Polytropic CEE" begin
+    for gamma in [1.4, 1.0, 5 / 3]
+        kappa = 0.5     # Scaling factor for the pressure.
+        equations = PolytropicEulerEquations2D(gamma, kappa)
+        u = SVector(1.1, -0.5, 2.34)
+
+        orientations = [1, 2]
+        for orientation in orientations
+            @test flux_lax_friedrichs(u, u, orientation, equations) ≈
+                  flux(u, orientation, equations)
+        end
+
+        normal_directions = [SVector(1.0, 0.0),
+            SVector(0.0, 1.0),
+            SVector(0.5, -0.5),
+            SVector(-1.2, 0.3)]
+
+        for normal_direction in normal_directions
+            @test flux_lax_friedrichs(u, u, normal_direction, equations) ≈
+                  flux(u, normal_direction, equations)
+        end
+    end
+end
+
 @timed_testset "Consistency check for HLL flux with Davis wave speed estimates: LEE" begin
     flux_hll = FluxHLL(min_max_speed_davis)