diff --git a/docs/literate/make.jl b/docs/literate/make.jl
index b620f85c975..a04d8a0b333 100644
--- a/docs/literate/make.jl
+++ b/docs/literate/make.jl
@@ -51,12 +51,25 @@ function create_tutorials(files)
     # Run tests on all tutorial files
     @testset "TrixiTutorials" begin
         for (i, (title, filename)) in enumerate(files)
+            # Evaluate each tutorial in its own module to avoid leaking of
+            # function/variable names, polluting the namespace of later tutorials
+            # by stuff defined in earlier tutorials.
             if filename isa Vector # Several files of one topic
                 for j in eachindex(filename)
-                    @testset "$(filename[j][2][2])" begin include(joinpath(repo_src, filename[j][2][1], filename[j][2][2])) end
+                    mod = gensym(filename[j][2][2])
+                    @testset "$(filename[j][2][2])" begin
+                        @eval module $mod
+                            include(joinpath($repo_src, $(filename[j][2][1]), $(filename[j][2][2])))
+                        end
+                    end
                 end
             else # Single files
-                @testset "$title" begin include(joinpath(repo_src, filename)) end
+                mod = gensym(title)
+                @testset "$title" begin
+                    @eval module $mod
+                        include(joinpath($repo_src, $filename))
+                    end
+                end
             end
         end
     end
diff --git a/docs/literate/src/files/custom_semidiscretization.jl b/docs/literate/src/files/custom_semidiscretization.jl
new file mode 100644
index 00000000000..fd432fb0826
--- /dev/null
+++ b/docs/literate/src/files/custom_semidiscretization.jl
@@ -0,0 +1,324 @@
+#src # Custom semidiscretizations
+
+# As described in the [overview section](@ref overview-semidiscretizations),
+# semidiscretizations are high-level descriptions of spatial discretizations
+# in Trixi.jl. Trixi.jl's main focus is on hyperbolic conservation
+# laws represented in a [`SemidiscretizationHyperbolic`](@ref).
+# Hyperbolic-parabolic problems based on the advection-diffusion equation or
+# the compressible Navier-Stokes equations can be represented in a
+# [`SemidiscretizationHyperbolicParabolic`](@ref). This is described in the
+# [basic tutorial on parabolic terms](@ref parabolic_terms) and its extension to
+# [custom parabolic terms](@ref adding_new_parabolic_terms).
+# In this tutorial, we will describe how these semidiscretizations work and how
+# they can be used to create custom semidiscretizations involving also other tasks.
+
+
+# ## Overview of the right-hand side evaluation
+
+# The semidiscretizations provided by Trixi.jl are set up to create `ODEProblem`s from the
+# [SciML ecosystem for ordinary differential equations](https://diffeq.sciml.ai/latest/).
+# In particular, a spatial semidiscretization can be wrapped in an ODE problem
+# using [`semidiscretize`](@ref), which returns an `ODEProblem`. This `ODEProblem`
+# bundles an initial condition, a right-hand side (RHS) function, the time span,
+# and possible parameters. The `ODEProblem`s created by Trixi.jl use the semidiscretization
+# passed to [`semidiscretize`](@ref) as a parameter.
+# For a [`SemidiscretizationHyperbolic`](@ref), the `ODEProblem` wraps
+# `Trixi.rhs!` as ODE RHS.
+# For a [`SemidiscretizationHyperbolicParabolic`](@ref),  Trixi.jl
+# uses a `SplitODEProblem` combining `Trixi.rhs_parabolic!` for the
+# (potentially) stiff part and `Trixi.rhs!` for the other part.
+
+
+# ## Standard Trixi.jl setup
+
+# In this tutorial, we will consider the linear advection equation
+# with source term
+# ```math
+# \partial_t u(t,x) + \partial_x u(t,x) = -\exp(-t) \sin\bigl(\pi (x - t) \bigr)
+# ```
+# with periodic boundary conditions in the domain `[-1, 1]` as a
+# model problem.
+# The initial condition is
+# ```math
+# u(0,x) = \sin(\pi x).
+# ```
+# The source term results in some damping and the analytical solution
+# ```math
+# u(t,x) = \exp(-t) \sin\bigl(\pi (x - t) \bigr).
+# ```
+# First, we discretize this equation using the standard functionality
+# of Trixi.jl.
+
+using Trixi, OrdinaryDiffEq, Plots
+
+# The linear scalar advection equation is already implemented in
+# Trixi.jl as [`LinearScalarAdvectionEquation1D`](@ref). We construct
+# it with an advection velocity `1.0`.
+
+equations = LinearScalarAdvectionEquation1D(1.0)
+
+# Next, we use a standard [`DGSEM`](@ref) solver.
+
+solver = DGSEM(polydeg = 3)
+
+# We create a simple [`TreeMesh`](@ref) in 1D.
+
+coordinates_min = (-1.0,)
+coordinates_max = (+1.0,)
+mesh = TreeMesh(coordinates_min, coordinates_max;
+                initial_refinement_level = 4,
+                n_cells_max = 10^4,
+                periodicity = true)
+
+# We wrap everything in in a semidiscretization and pass the source
+# terms as a standard Julia function. Please note that Trixi.jl uses
+# `SVector`s from
+# [StaticArrays.jl](https://github.com/JuliaArrays/StaticArrays.jl)
+# to store the conserved variables `u`. Thus, the return value of the
+# source terms must be wrapped in an `SVector` - even if we consider
+# just a scalar problem.
+
+function initial_condition(x, t, equations)
+    return SVector(exp(-t) * sinpi(x[1] - t))
+end
+
+function source_terms_standard(u, x, t, equations)
+    return -initial_condition(x, t, equations)
+end
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition,
+                                    solver;
+                                    source_terms = source_terms_standard)
+
+# Now, we can create the `ODEProblem`, solve the resulting ODE
+# using a time integration method from
+# [OrdinaryDiffEq.jl](https://github.com/SciML/OrdinaryDiffEq.jl),
+# and visualize the numerical solution at the final time using
+# [Plots.jl](https://github.com/JuliaPlots/Plots.jl).
+
+tspan = (0.0, 3.0)
+ode = semidiscretize(semi, tspan)
+
+sol = solve(ode, RDPK3SpFSAL49(); ode_default_options()...)
+
+plot(sol; label = "numerical sol.", legend = :topright)
+
+# We can also plot the analytical solution for comparison.
+# Since Trixi.jl uses `SVector`s for the variables, we take their `first`
+# (and only) component to get the scalar value for manual plotting.
+
+let
+   x = range(-1.0, 1.0; length = 200)
+   plot!(x, first.(initial_condition.(x, sol.t[end], equations)),
+         label = "analytical sol.", linestyle = :dash, legend = :topright)
+end
+
+# We can also add the initial condition to the plot.
+
+plot!(sol.u[1], semi, label = "u0", linestyle = :dot, legend = :topleft)
+
+# You can of course also use some
+# [callbacks](https://trixi-framework.github.io/Trixi.jl/stable/callbacks/)
+# provided by Trixi.jl as usual.
+
+summary_callback = SummaryCallback()
+analysis_interval = 100
+analysis_callback = AnalysisCallback(semi; interval = analysis_interval)
+alive_callback = AliveCallback(; analysis_interval)
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback,
+                        alive_callback)
+
+sol = solve(ode, RDPK3SpFSAL49();
+            ode_default_options()..., callback = callbacks)
+summary_callback()
+
+
+# ## Using a custom ODE right-hand side function
+
+# Next, we will solve the same problem but use our own ODE RHS function.
+# To demonstrate this, we will artificially create a global variable
+# containing the current time of the simulation.
+
+const GLOBAL_TIME = Ref(0.0)
+
+function source_terms_custom(u, x, t, equations)
+    t = GLOBAL_TIME[]
+    return -initial_condition(x, t, equations)
+end
+
+# Next, we create our own RHS function to update the global time of
+# the simulation before calling the RHS function from Trixi.jl.
+
+function rhs_source_custom!(du_ode, u_ode, semi, t)
+    GLOBAL_TIME[] = t
+    Trixi.rhs!(du_ode, u_ode, semi, t)
+end
+
+# Next, we create an `ODEProblem` manually copying over the data from
+# the one we got from [`semidiscretize`](@ref) earlier.
+
+ode_source_custom = ODEProblem(rhs_source_custom!,
+                               ode.u0,
+                               ode.tspan,
+                               ode.p #= semi =#)
+sol_source_custom = solve(ode_source_custom, RDPK3SpFSAL49();
+                          ode_default_options()...)
+
+plot(sol_source_custom; label = "numerical sol.")
+let
+    x = range(-1.0, 1.0; length = 200)
+    plot!(x, first.(initial_condition.(x, sol_source_custom.t[end], equations)),
+          label = "analytical sol.", linestyle = :dash, legend = :topleft)
+end
+plot!(sol_source_custom.u[1], semi, label = "u0", linestyle = :dot, legend = :topleft)
+
+# This also works with callbacks as usual.
+
+summary_callback = SummaryCallback()
+analysis_interval = 100
+analysis_callback = AnalysisCallback(semi; interval = analysis_interval)
+alive_callback = AliveCallback(; analysis_interval)
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback,
+                        alive_callback)
+
+sol = solve(ode_source_custom, RDPK3SpFSAL49();
+            ode_default_options()..., callback = callbacks)
+summary_callback()
+
+
+# ## Setting up a custom semidiscretization
+
+# Using a global constant is of course not really nice from a software
+# engineering point of view. Thus, it can often be useful to collect
+# additional data in the parameters of the `ODEProblem`. Thus, it is
+# time to create our own semidiscretization. Here, we create a small
+# wrapper of a standard semidiscretization of Trixi.jl and the current
+# global time of the simulation.
+
+struct CustomSemidiscretization{Semi, T} <: Trixi.AbstractSemidiscretization
+    semi::Semi
+    t::T
+end
+
+semi_custom = CustomSemidiscretization(semi, Ref(0.0))
+
+# To get pretty printing in the REPL, you can consider specializing
+#
+# - `Base.show(io::IO, parameters::CustomSemidiscretization)`
+# - `Base.show(io::IO, ::MIME"text/plain", parameters::CustomSemidiscretization)`
+#
+# for your custom semidiscretiation.
+
+# Next, we create our own source terms that use the global time stored
+# in the custom semidiscretiation.
+
+source_terms_custom_semi = let semi_custom = semi_custom
+    function source_terms_custom_semi(u, x, t, equations)
+        t = semi_custom.t[]
+        return -initial_condition(x, t, equations)
+    end
+end
+
+# We also create a custom ODE RHS to update the current global time
+# stored in the custom semidiscretization. We unpack the standard
+# semidiscretization created by Trixi.jl and pass it to `Trixi.rhs!`.
+
+function rhs_semi_custom!(du_ode, u_ode, semi_custom, t)
+    semi_custom.t[] = t
+    Trixi.rhs!(du_ode, u_ode, semi_custom.semi, t)
+end
+
+# Finally, we set up an `ODEProblem` and solve it numerically.
+
+ode_semi_custom = ODEProblem(rhs_semi_custom!,
+                             ode.u0,
+                             ode.tspan,
+                             semi_custom)
+sol_semi_custom = solve(ode_semi_custom, RDPK3SpFSAL49();
+                        ode_default_options()...)
+
+# If we want to make use of additional functionality provided by
+# Trixi.jl, e.g., for plotting, we need to implement a few additional
+# specializations. In this case, we forward everything to the standard
+# semidiscretization provided by Trixi.jl wrapped in our custom
+# semidiscretization.
+
+Base.ndims(semi::CustomSemidiscretization) = ndims(semi.semi)
+function Trixi.mesh_equations_solver_cache(semi::CustomSemidiscretization)
+    Trixi.mesh_equations_solver_cache(semi.semi)
+end
+
+# Now, we can plot the numerical solution as usual.
+
+plot(sol_semi_custom; label = "numerical sol.")
+let
+    x = range(-1.0, 1.0; length = 200)
+    plot!(x, first.(initial_condition.(x, sol_semi_custom.t[end], equations)),
+          label = "analytical sol.", linestyle = :dash, legend = :topleft)
+end
+plot!(sol_semi_custom.u[1], semi, label = "u0", linestyle = :dot, legend = :topleft)
+
+# This also works with many callbacks as usual. However, the
+# [`AnalysisCallback`](@ref) requires some special handling since it
+# makes use of a performance counter contained in the standard
+# semidiscretizations of Trixi.jl to report some
+# [performance metrics](@ref performance-metrics).
+# Here, we forward all accesses to the performance counter to the
+# wrapped semidiscretization.
+
+function Base.getproperty(semi::CustomSemidiscretization, s::Symbol)
+    if s === :performance_counter
+        wrapped_semi = getfield(semi, :semi)
+        wrapped_semi.performance_counter
+    else
+        getfield(semi, s)
+    end
+end
+
+# Moreover, the [`AnalysisCallback`](@ref) also performs some error
+# calculations. We also need to forward them to the wrapped
+# semidiscretization.
+
+function Trixi.calc_error_norms(func, u, t, analyzer,
+                                semi::CustomSemidiscretization,
+                                cache_analysis)
+    Trixi.calc_error_norms(func, u, t, analyzer,
+                           semi.semi,
+                           cache_analysis)
+end
+
+# Now, we can work with the callbacks used before as usual.
+
+summary_callback = SummaryCallback()
+analysis_interval = 100
+analysis_callback = AnalysisCallback(semi_custom;
+                                     interval = analysis_interval)
+alive_callback = AliveCallback(; analysis_interval)
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback,
+                        alive_callback)
+
+sol = solve(ode_semi_custom, RDPK3SpFSAL49();
+            ode_default_options()..., callback = callbacks)
+summary_callback()
+
+# For even more advanced usage of custom semidiscretizations, you
+# may look at the source code of the ones contained in Trixi.jl, e.g.,
+# - [`SemidiscretizationHyperbolicParabolic`](@ref)
+# - [`SemidiscretizationEulerGravity`](@ref)
+# - [`SemidiscretizationEulerAcoustics`](@ref)
+# - [`SemidiscretizationCoupled`](@ref)
+
+
+# ## Package versions
+
+# These results were obtained using the following versions.
+
+using InteractiveUtils
+versioninfo()
+
+using Pkg
+Pkg.status(["Trixi", "OrdinaryDiffEq", "Plots"],
+           mode=PKGMODE_MANIFEST)
diff --git a/docs/literate/src/files/index.jl b/docs/literate/src/files/index.jl
index 0c8de66bf42..d42695611f6 100644
--- a/docs/literate/src/files/index.jl
+++ b/docs/literate/src/files/index.jl
@@ -76,21 +76,30 @@
 # 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.
 
-# ### [10 Adaptive mesh refinement](@ref adaptive_mesh_refinement)
+# ### [10 Parabolic terms](@ref parabolic_terms)
+#-
+# This tutorial describes how parabolic terms are implemented in Trixi.jl, e.g.,
+# to solve the advection-diffusion equation.
+
+# ### [11 Adding new parabolic terms](@ref adding_new_parabolic_terms)
+#-
+# This tutorial describes how new parabolic terms can be implemented using Trixi.jl.
+
+# ### [12 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.
 
-# ### [11 Structured mesh with curvilinear mapping](@ref structured_mesh_mapping)
+# ### [13 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.
 
-# ### [12 Unstructured meshes with HOHQMesh.jl](@ref hohqmesh_tutorial)
+# ### [14 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
@@ -99,19 +108,24 @@
 # 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`.
 
-# ### [13 Explicit time stepping](@ref time_stepping)
+# ### [15 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.
 
-# ### [14 Differentiable programming](@ref differentiable_programming)
+# ### [16 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.
 
+# ### [17 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.
+
 
 
 # ## Examples in Trixi.jl
diff --git a/docs/make.jl b/docs/make.jl
index 57629577ddb..f882fcf1219 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -68,6 +68,7 @@ files = [
     # Topic: other stuff
     "Explicit time stepping" => "time_stepping.jl",
     "Differentiable programming" => "differentiable_programming.jl",
+    "Custom semidiscretizations" => "custom_semidiscretization.jl"
     ]
 tutorials = create_tutorials(files)
 
diff --git a/docs/src/callbacks.md b/docs/src/callbacks.md
index 1d3e5e34b51..7f44dfd5925 100644
--- a/docs/src/callbacks.md
+++ b/docs/src/callbacks.md
@@ -30,7 +30,7 @@ An example elixir using AMR can be found at [`examples/tree_2d_dgsem/elixir_adve
 The [`AnalysisCallback`](@ref) can be used to analyze the numerical solution, e.g. calculate
 errors or user-specified integrals, and print the results to the screen. The results can also be
 saved in a file. An example can be found at [`examples/tree_2d_dgsem/elixir_euler_vortex.jl`](https://github.com/trixi-framework/Trixi.jl/blob/main/examples/tree_2d_dgsem/elixir_euler_vortex.jl).
-In [Performance metrics of the `AnalysisCallback`](@ref) you can find a detailed
+In [Performance metrics of the `AnalysisCallback`](@ref performance-metrics) you can find a detailed
 description of the different performance metrics the `AnalysisCallback` computes.
 
 ### I/O
@@ -106,7 +106,7 @@ will yield the following plot:
   the automated performance measurements, including an output of the recorded timers after a simulation.
 * The [`VisualizationCallback`](@ref) can be used for in-situ visualization. See
   [Visualizing results during a simulation](@ref).
-* The [`TrivialCallback`](@ref) does nothing and can be used to to easily disable some callbacks
+* The [`TrivialCallback`](@ref) does nothing and can be used to easily disable some callbacks
   via [`trixi_include`](@ref).
 
 ### Equation-specific callbacks
diff --git a/docs/src/overview.md b/docs/src/overview.md
index 46bc28b6025..51a6272ae8e 100644
--- a/docs/src/overview.md
+++ b/docs/src/overview.md
@@ -16,7 +16,7 @@ to solve a PDE numerically are the spatial semidiscretization and the time
 integration scheme.
 
 
-## Semidiscretizations
+## [Semidiscretizations](@id overview-semidiscretizations)
 
 Semidiscretizations are high-level descriptions of spatial discretizations
 specialized for certain PDEs. Trixi.jl's main focus is on hyperbolic conservation
diff --git a/docs/src/parallelization.md b/docs/src/parallelization.md
index d56777c9af4..e55471bb256 100644
--- a/docs/src/parallelization.md
+++ b/docs/src/parallelization.md
@@ -22,7 +22,7 @@ julia --threads=4
 If both the environment variable and the command line argument are specified at
 the same time, the latter takes precedence.
 
-If you use time integration methods from 
+If you use time integration methods from
 [OrdinaryDiffEq.jl](https://github.com/SciML/OrdinaryDiffEq.jl)
 and want to use multiple threads therein, you need to set the keyword argument
 `thread=OrdinaryDiffEq.True()` of the algorithms, as described in the
@@ -143,7 +143,7 @@ To start Trixi.jl in parallel with MPI, there are three options:
    Switching between panes can be done by `Ctrl+b` followed by `o`.
    As of March 2022, newer versions of tmpi also support mpich, which is the default
    backend of MPI.jl (via MPICH_Jll.jl). To use this setup, you need to install
-   `mpiexecjl` as described in the 
+   `mpiexecjl` as described in the
    [documentation of MPI.jl](https://juliaparallel.org/MPI.jl/v0.20/usage/#Julia-wrapper-for-mpiexec)
    and make it available as `mpirun`, e.g., via a symlink of the form
    ```bash
@@ -161,7 +161,7 @@ To start Trixi.jl in parallel with MPI, there are three options:
 
 ### [Performance](@id parallel_performance)
 For information on how to evaluate the parallel performance of Trixi.jl, please
-have a look at the [Performance metrics of the `AnalysisCallback`](@ref)
+have a look at the [Performance metrics of the `AnalysisCallback`](@ref performance-metrics)
 section, specifically at the descriptions of the performance index (PID).
 
 
diff --git a/docs/src/performance.md b/docs/src/performance.md
index 428672ec75f..bbe3a3390b7 100644
--- a/docs/src/performance.md
+++ b/docs/src/performance.md
@@ -170,7 +170,7 @@ As a rule of thumb:
 - Consider using `@nospecialize` for methods like custom implementations of `Base.show`.
 
 
-## Performance metrics of the `AnalysisCallback`
+## [Performance metrics of the `AnalysisCallback`](@id performance-metrics)
 The [`AnalysisCallback`](@ref) computes two performance indicators that you can use to
 evaluate the serial and parallel performance of Trixi.jl. They represent
 measured run times that are normalized by the number of `rhs!` evaluations and