diff --git a/NEWS.md b/NEWS.md
index 873b8c1c162..6417045e45f 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -37,5 +37,6 @@ for human readability.
   `cons2prim`
 - `varnames_cons(equations)` → `varnames(cons2cons, equations)`
 - `varnames_prim(equations)` → `varnames(cons2prim, equations)`
+- `max_dt(u, t, mesh, ...)` → `max_dt_hyperbolic(u, t, mesh, ...)` ([#378](https://github.com/trixi-framework/Trixi.jl/pull/378))
 
 #### Removed
diff --git a/examples/2d/elixir_advdiff_amr.jl b/examples/2d/elixir_advdiff_amr.jl
new file mode 100644
index 00000000000..ce0a4907e38
--- /dev/null
+++ b/examples/2d/elixir_advdiff_amr.jl
@@ -0,0 +1,74 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the linear advection equation
+
+advectionvelocity = (1.8, 0.6)
+nu = 2.4e-1
+equations = LinearAdvectionDiffusionEquation2D(advectionvelocity, nu)
+
+# Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux
+polydeg = 3
+solver_hyperbolic = DGSEM(polydeg, flux_lax_friedrichs)
+
+coordinates_min = (-4, -4)
+coordinates_max = ( 4,  4)
+
+# Create a uniformely refined mesh with periodic boundaries
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level=3,
+                n_cells_max=30_000) # set maximum capacity of tree data structure
+
+# A semidiscretization collects data structures and functions for the spatial discretization
+initial_condition = initial_condition_gauss
+semi = SemidiscretizationHyperbolicParabolicBR1(mesh, equations, initial_condition, solver_hyperbolic)
+
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+# Create ODE problem with time span from 0.0 to 1.0
+tspan = (0.0, 1.0)
+ode = semidiscretize(semi, tspan);
+
+# At the beginning of the main loop, the SummaryCallback prints a summary of the simulation setup
+# and resets the timers
+summary_callback = SummaryCallback()
+
+# The AnalysisCallback allows to analyse the solution in regular intervals and prints the results
+analysis_callback = AnalysisCallback(semi, interval=100)
+
+# The SaveSolutionCallback allows to save the solution to a file in regular intervals
+save_solution = SaveSolutionCallback(interval=100,
+                                     solution_variables=cons2prim)
+
+amr_controller = ControllerThreeLevel(semi, IndicatorMax(semi, variable=first),
+                                      base_level=3,
+                                      med_level=4, med_threshold=0.1,
+                                      max_level=5, max_threshold=0.5)
+amr_callback = AMRCallback(semi, amr_controller,
+                           interval=5,
+                           adapt_initial_condition=true,
+                           adapt_initial_condition_only_refine=true)
+
+# The StepsizeCallback handles the re-calculcation of the maximum Δt after each time step
+# The first CFL number is for the hyperbolic system, the second for the parabolic system
+stepsize_callback = StepsizeCallback(cfl=(1.6, 0.5))
+
+# Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver
+callbacks = CallbackSet(summary_callback, analysis_callback, save_solution, amr_callback,
+                        stepsize_callback)
+
+
+###############################################################################
+# run the simulation
+
+# OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
+            dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
+            save_everystep=false, callback=callbacks);
+
+# Print the timer summary
+summary_callback()
diff --git a/examples/2d/elixir_advdiff_basic.jl b/examples/2d/elixir_advdiff_basic.jl
new file mode 100644
index 00000000000..ac5b44f648c
--- /dev/null
+++ b/examples/2d/elixir_advdiff_basic.jl
@@ -0,0 +1,68 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the linear advection equation
+
+advectionvelocity = (0.2, 0.6)
+nu = 1.2e-2
+equations = LinearAdvectionDiffusionEquation2D(advectionvelocity, nu)
+
+# Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux
+polydeg = 4
+solver_hyperbolic = DGSEM(polydeg, flux_lax_friedrichs)
+
+coordinates_min = (-1, -1) # minimum coordinates (min(x), min(y))
+coordinates_max = ( 1,  1) # maximum coordinates (max(x), max(y))
+refinement_patches = (
+  (type="box", coordinates_min=(0.0, -1.0), coordinates_max=(1.0, 1.0)),
+)
+
+# Create a uniformely refined mesh with periodic boundaries
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level=2,
+                refinement_patches=refinement_patches,
+                n_cells_max=30_000) # set maximum capacity of tree data structure
+
+# A semidiscretization collects data structures and functions for the spatial discretization
+initial_condition = initial_condition_convergence_test
+semi = SemidiscretizationHyperbolicParabolicBR1(mesh, equations, initial_condition, solver_hyperbolic)
+
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+# Create ODE problem with time span from 0.0 to 1.0
+tspan = (0.0, 1.0)
+ode = semidiscretize(semi, tspan);
+
+# At the beginning of the main loop, the SummaryCallback prints a summary of the simulation setup
+# and resets the timers
+summary_callback = SummaryCallback()
+
+# The AnalysisCallback allows to analyse the solution in regular intervals and prints the results
+analysis_callback = AnalysisCallback(semi, interval=100)
+
+# The SaveSolutionCallback allows to save the solution to a file in regular intervals
+save_solution = SaveSolutionCallback(interval=100,
+                                     solution_variables=cons2prim)
+
+# The StepsizeCallback handles the re-calculcation of the maximum Δt after each time step
+# The first CFL number is for the hyperbolic system, the second for the parabolic system
+stepsize_callback = StepsizeCallback(cfl=(1.0, 0.1))
+
+# Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver
+callbacks = CallbackSet(summary_callback, analysis_callback, save_solution, stepsize_callback)
+
+
+###############################################################################
+# run the simulation
+
+# OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
+            dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
+            save_everystep=false, callback=callbacks);
+
+# Print the timer summary
+summary_callback()
diff --git a/examples/2d/elixir_heat_amr.jl b/examples/2d/elixir_heat_amr.jl
new file mode 100644
index 00000000000..00981bafbb3
--- /dev/null
+++ b/examples/2d/elixir_heat_amr.jl
@@ -0,0 +1,71 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the linear advection equation
+
+nu = 4.8e-1
+equations = HeatEquation2D(nu)
+
+# Create DG solver with polynomial degree = 3 and the central flux as surface flux
+solver = DGSEM(3, flux_central)
+
+coordinates_min = (-4, -4) # minimum coordinates (min(x), min(y))
+coordinates_max = ( 4,  4) # maximum coordinates (max(x), max(y))
+
+# Create a uniformely refined mesh with periodic boundaries
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level=3,
+                n_cells_max=30_000) # set maximum capacity of tree data structure
+
+# A semidiscretization collects data structures and functions for the spatial discretization
+initial_condition = initial_condition_gauss
+semi = SemidiscretizationParabolicAuxVars(mesh, equations, initial_condition, solver)
+
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+# Create ODE problem with time span from 0.0 to 1.0
+tspan = (0.0, 1.0)
+ode = semidiscretize(semi, tspan);
+
+# At the beginning of the main loop, the SummaryCallback prints a summary of the simulation setup
+# and resets the timers
+summary_callback = SummaryCallback()
+
+# The AnalysisCallback allows to analyse the solution in regular intervals and prints the results
+analysis_callback = AnalysisCallback(semi, interval=100)
+
+# The SaveSolutionCallback allows to save the solution to a file in regular intervals
+save_solution = SaveSolutionCallback(interval=100,
+                                     solution_variables=cons2prim)
+
+amr_controller = ControllerThreeLevel(semi, IndicatorMax(semi, variable=first),
+                                      base_level=3,
+                                      med_level=4, med_threshold=0.05,
+                                      max_level=5, max_threshold=0.3)
+amr_callback = AMRCallback(semi, amr_controller,
+                           interval=5,
+                           adapt_initial_condition=true,
+                           adapt_initial_condition_only_refine=true)
+
+# The StepsizeCallback handles the re-calculcation of the maximum Δt after each time step
+stepsize_callback = StepsizeCallback(cfl=0.5)
+
+# Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver
+callbacks = CallbackSet(summary_callback, analysis_callback, save_solution, amr_callback,
+                        stepsize_callback)
+
+
+###############################################################################
+# run the simulation
+
+# OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
+            dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
+            save_everystep=false, callback=callbacks, maxiters=1e5);
+
+# Print the timer summary
+summary_callback()
diff --git a/examples/2d/elixir_heat_basic.jl b/examples/2d/elixir_heat_basic.jl
new file mode 100644
index 00000000000..18d0be4cf01
--- /dev/null
+++ b/examples/2d/elixir_heat_basic.jl
@@ -0,0 +1,61 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the linear advection equation
+
+nu = 1.2e-2
+equations = HeatEquation2D(nu)
+
+# Create DG solver with polynomial degree = 3 and the central flux as surface flux
+solver = DGSEM(3, flux_central)
+
+coordinates_min = (-1, -1) # minimum coordinates (min(x), min(y))
+coordinates_max = ( 1,  1) # maximum coordinates (max(x), max(y))
+
+# Create a uniformely refined mesh with periodic boundaries
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level=2,
+                n_cells_max=30_000) # set maximum capacity of tree data structure
+
+# A semidiscretization collects data structures and functions for the spatial discretization
+initial_condition = initial_condition_convergence_test
+semi = SemidiscretizationParabolicAuxVars(mesh, equations, initial_condition, solver)
+
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+# Create ODE problem with time span from 0.0 to 50.0
+tspan = (0.0, 50.0)
+ode = semidiscretize(semi, tspan);
+
+# At the beginning of the main loop, the SummaryCallback prints a summary of the simulation setup
+# and resets the timers
+summary_callback = SummaryCallback()
+
+# The AnalysisCallback allows to analyse the solution in regular intervals and prints the results
+analysis_callback = AnalysisCallback(semi, interval=100)
+
+# The SaveSolutionCallback allows to save the solution to a file in regular intervals
+save_solution = SaveSolutionCallback(interval=100,
+                                     solution_variables=cons2prim)
+
+# The StepsizeCallback handles the re-calculcation of the maximum Δt after each time step
+stepsize_callback = StepsizeCallback(cfl=0.5)
+
+# Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver
+callbacks = CallbackSet(summary_callback, analysis_callback, save_solution, stepsize_callback)
+
+
+###############################################################################
+# run the simulation
+
+# OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
+            dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
+            save_everystep=false, callback=callbacks, maxiters=1e5);
+
+# Print the timer summary
+summary_callback()
diff --git a/examples/2d/elixir_heat_nonperiodic.jl b/examples/2d/elixir_heat_nonperiodic.jl
new file mode 100644
index 00000000000..d6c63b90bca
--- /dev/null
+++ b/examples/2d/elixir_heat_nonperiodic.jl
@@ -0,0 +1,69 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the linear advection equation
+
+nu = 1.2e-2
+equations = HeatEquation2D(nu)
+
+# Create DG solver with polynomial degree = 3 and the central flux as surface flux
+solver = DGSEM(3, flux_central)
+
+coordinates_min = (-1, -1) # minimum coordinates (min(x), min(y))
+coordinates_max = ( 1,  1) # maximum coordinates (max(x), max(y))
+
+# Create a uniformely refined mesh with periodic boundaries
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level=2,
+                n_cells_max=30_000, # set maximum capacity of tree data structure
+                periodicity=(false, true))
+
+# A semidiscretization collects data structures and functions for the spatial discretization
+initial_condition = initial_condition_sin_x
+
+# 1 => -x, 2 => +x, 3 => -y, 4 => +y as usual for orientations
+boundary_conditions = (x_neg=boundary_condition_sin_x,
+                       x_pos=boundary_condition_sin_x,
+                       y_neg=boundary_condition_periodic,
+                       y_pos=boundary_condition_periodic)
+semi = SemidiscretizationParabolicAuxVars(mesh, equations, initial_condition, solver,
+                                          boundary_conditions=boundary_conditions)
+
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+# Create ODE problem with time span from 0.0 to 100.0
+tspan = (0.0, 100.0)
+ode = semidiscretize(semi, tspan);
+
+# At the beginning of the main loop, the SummaryCallback prints a summary of the simulation setup
+# and resets the timers
+summary_callback = SummaryCallback()
+
+# The AnalysisCallback allows to analyse the solution in regular intervals and prints the results
+analysis_callback = AnalysisCallback(semi, interval=100)
+
+# The SaveSolutionCallback allows to save the solution to a file in regular intervals
+save_solution = SaveSolutionCallback(interval=100,
+                                     solution_variables=cons2prim)
+
+# The StepsizeCallback handles the re-calculcation of the maximum Δt after each time step
+stepsize_callback = StepsizeCallback(cfl=0.5)
+
+# Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver
+callbacks = CallbackSet(summary_callback, analysis_callback, save_solution, stepsize_callback)
+
+
+###############################################################################
+# run the simulation
+
+# OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
+            dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
+            save_everystep=false, callback=callbacks, maxiters=1e5);
+
+# Print the timer summary
+summary_callback()
diff --git a/examples/2d/elixir_heat_source_terms.jl b/examples/2d/elixir_heat_source_terms.jl
new file mode 100644
index 00000000000..e3b3a87944b
--- /dev/null
+++ b/examples/2d/elixir_heat_source_terms.jl
@@ -0,0 +1,62 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the linear advection equation
+
+nu = 1.2e-2
+equations = HeatEquation2D(nu)
+
+# Create DG solver with polynomial degree = 3 and the central flux as surface flux
+solver = DGSEM(3, flux_central)
+
+coordinates_min = (-1, -1) # minimum coordinates (min(x), min(y))
+coordinates_max = ( 1,  1) # maximum coordinates (max(x), max(y))
+
+# Create a uniformely refined mesh with periodic boundaries
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level=2,
+                n_cells_max=30_000) # set maximum capacity of tree data structure
+
+# A semidiscretization collects data structures and functions for the spatial discretization
+initial_condition = initial_condition_poisson_periodic
+semi = SemidiscretizationParabolicAuxVars(mesh, equations, initial_condition, solver,
+                                          source_terms=source_terms_poisson_periodic)
+
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+# Create ODE problem with time span from 0.0 to 100.0
+tspan = (0.0, 100.0)
+ode = semidiscretize(semi, tspan);
+
+# At the beginning of the main loop, the SummaryCallback prints a summary of the simulation setup
+# and resets the timers
+summary_callback = SummaryCallback()
+
+# The AnalysisCallback allows to analyse the solution in regular intervals and prints the results
+analysis_callback = AnalysisCallback(semi, interval=100)
+
+# The SaveSolutionCallback allows to save the solution to a file in regular intervals
+save_solution = SaveSolutionCallback(interval=100,
+                                     solution_variables=cons2prim)
+
+# The StepsizeCallback handles the re-calculcation of the maximum Δt after each time step
+stepsize_callback = StepsizeCallback(cfl=0.5)
+
+# Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver
+callbacks = CallbackSet(summary_callback, analysis_callback, save_solution, stepsize_callback)
+
+
+###############################################################################
+# run the simulation
+
+# OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
+            dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
+            save_everystep=false, callback=callbacks, maxiters=1e5);
+
+# Print the timer summary
+summary_callback()
diff --git a/src/Trixi.jl b/src/Trixi.jl
index 332be6a7922..24882d43298 100644
--- a/src/Trixi.jl
+++ b/src/Trixi.jl
@@ -64,6 +64,8 @@ include("semidiscretization/semidiscretization_hyperbolic.jl")
 include("callbacks_step/callbacks_step.jl")
 include("callbacks_stage/callbacks_stage.jl")
 include("semidiscretization/semidiscretization_euler_gravity.jl")
+include("semidiscretization/semidiscretization_parabolic_auxvars.jl")
+include("semidiscretization/semidiscretization_hyperbolic_parabolic.jl")
 include("time_integration/time_integration.jl")
 
 # `trixi_include` and special elixirs such as `convergence_test`
@@ -83,7 +85,9 @@ export AcousticPerturbationEquations2D,
        HyperbolicDiffusionEquations1D, HyperbolicDiffusionEquations2D, HyperbolicDiffusionEquations3D,
        LinearScalarAdvectionEquation1D, LinearScalarAdvectionEquation2D, LinearScalarAdvectionEquation3D,
        InviscidBurgersEquation1D,
-       LatticeBoltzmannEquations2D, LatticeBoltzmannEquations3D
+       LatticeBoltzmannEquations2D, LatticeBoltzmannEquations3D,
+       HeatEquation2D,
+       LinearAdvectionDiffusionEquation2D
 
 export flux, flux_central, flux_lax_friedrichs, flux_hll, flux_hllc, flux_godunov,
        flux_chandrashekar, flux_ranocha, flux_derigs_etal, flux_kennedy_gruber, flux_shima_etal,
@@ -130,6 +134,7 @@ export initial_condition_lid_driven_cavity, boundary_condition_lid_driven_cavity
 export initial_condition_couette_steady, initial_condition_couette_unsteady, boundary_condition_couette
 export initial_condition_gauss_wall
 export initial_condition_monopole, boundary_condition_monopole
+export initial_condition_sin_x, boundary_condition_sin_x
 
 export cons2cons, cons2prim, prim2cons, cons2macroscopic, cons2state, cons2mean,
        cons2entropy, entropy2cons
@@ -148,7 +153,9 @@ export DG,
 export nelements, nnodes, nvariables,
        eachelement, eachnode, eachvariable
 
-export SemidiscretizationHyperbolic, semidiscretize, compute_coefficients, integrate
+export SemidiscretizationHyperbolic, SemidiscretizationParabolicAuxVars,
+       SemidiscretizationHyperbolicParabolic, SemidiscretizationHyperbolicParabolicBR1,
+       semidiscretize, compute_coefficients, integrate
 
 export SemidiscretizationEulerGravity, ParametersEulerGravity,
        timestep_gravity_erk52_3Sstar!, timestep_gravity_carpenter_kennedy_erk54_2N!
diff --git a/src/callbacks_step/stepsize.jl b/src/callbacks_step/stepsize.jl
index d2f284454b3..88f697ea82a 100644
--- a/src/callbacks_step/stepsize.jl
+++ b/src/callbacks_step/stepsize.jl
@@ -1,12 +1,13 @@
 
 """
-    StepsizeCallback(; cfl=1.0)
+    StepsizeCallback(; cfl)
 
 Set the time step size according to a CFL condition with CFL number `cfl`
-if the time integration method isn't adaptive itself.
+if the time integration method isn't adaptive itself. When using a semidiscretization that wraps
+multiple solvers, you might need to provide a tuple of CFL numbers.
 """
-mutable struct StepsizeCallback{RealT}
-  cfl_number::RealT
+mutable struct StepsizeCallback{CFLNumber}
+  cfl_number::CFLNumber
 end
 
 
@@ -34,7 +35,7 @@ function Base.show(io::IO, ::MIME"text/plain", cb::DiscreteCallback{<:Any, <:Ste
 end
 
 
-function StepsizeCallback(; cfl::Real=1.0)
+function StepsizeCallback(; cfl)
 
   stepsize_callback = StepsizeCallback(cfl)
 
@@ -62,13 +63,10 @@ end
     t = integrator.t
     u_ode = integrator.u
     semi = integrator.p
-    mesh, equations, solver, cache = mesh_equations_solver_cache(semi)
     @unpack cfl_number = stepsize_callback
-    u = wrap_array(u_ode, mesh, equations, solver, cache)
 
-    dt = @timeit_debug timer() "calculate dt" cfl_number * max_dt(u, t, mesh,
-                                                                  have_constant_speed(equations), equations,
-                                                                  solver, cache)
+    dt = @timeit_debug timer() "calculate dt" max_dt(u_ode, t, cfl_number, semi)
+
     set_proposed_dt!(integrator, dt)
     integrator.opts.dtmax = dt
     integrator.dtcache = dt
@@ -90,13 +88,15 @@ function (cb::DiscreteCallback{Condition,Affect!})(ode::ODEProblem) where {Condi
   u_ode = ode.u0
   t = first(ode.tspan)
   semi = ode.p
-  mesh, equations, solver, cache = mesh_equations_solver_cache(semi)
-  u = wrap_array(u_ode, mesh, equations, solver, cache)
 
-  return cfl_number * max_dt(u, t, mesh, have_constant_speed(equations), equations, solver, cache)
+  return max_dt(u_ode, t, cfl_number, semi)
 end
 
 
+# FIXME: Deprecations introduced in v0.3
+@deprecate max_dt(u::AbstractArray, t, mesh::AbstractMesh, constant_speed, equations, dg::DG, cache) max_dt_hyperbolic(u, t, mesh, constant_speed, equations, dg, cache)
+
+
 include("stepsize_dg1d.jl")
 include("stepsize_dg2d.jl")
 include("stepsize_dg3d.jl")
diff --git a/src/callbacks_step/stepsize_dg1d.jl b/src/callbacks_step/stepsize_dg1d.jl
index 51a9cf5445e..56481d91159 100644
--- a/src/callbacks_step/stepsize_dg1d.jl
+++ b/src/callbacks_step/stepsize_dg1d.jl
@@ -1,6 +1,6 @@
 
-function max_dt(u::AbstractArray{<:Any,3}, t, mesh::Union{TreeMesh{1},StructuredMesh{1}},
-                constant_speed::Val{false}, equations, dg::DG, cache)
+function max_dt_hyperbolic(u::AbstractArray{<:Any,3}, t, mesh::Union{TreeMesh{1},StructuredMesh{1}},
+                           constant_speed::Val{false}, equations, dg::DG, cache)
   # to avoid a division by zero if the speed vanishes everywhere,
   # e.g. for steady-state linear advection
   max_scaled_speed = nextfloat(zero(t))
@@ -20,8 +20,8 @@ function max_dt(u::AbstractArray{<:Any,3}, t, mesh::Union{TreeMesh{1},Structured
 end
 
 
-function max_dt(u::AbstractArray{<:Any,3}, t, mesh::Union{TreeMesh{1},StructuredMesh{1}},
-                constant_speed::Val{true}, equations, dg::DG, cache)
+function max_dt_hyperbolic(u::AbstractArray{<:Any,3}, t, mesh::Union{TreeMesh{1},StructuredMesh{1}},
+                           constant_speed::Val{true}, equations, dg::DG, cache)
   # to avoid a division by zero if the speed vanishes everywhere,
   # e.g. for steady-state linear advection
   max_scaled_speed = nextfloat(zero(t))
diff --git a/src/callbacks_step/stepsize_dg2d.jl b/src/callbacks_step/stepsize_dg2d.jl
index aff86c74217..7604067966e 100644
--- a/src/callbacks_step/stepsize_dg2d.jl
+++ b/src/callbacks_step/stepsize_dg2d.jl
@@ -1,6 +1,6 @@
 
-function max_dt(u::AbstractArray{<:Any,4}, t, mesh::TreeMesh{2},
-                constant_speed::Val{false}, equations, dg::DG, cache)
+function max_dt_hyperbolic(u::AbstractArray{<:Any,4}, t, mesh::TreeMesh{2},
+                           constant_speed::Val{false}, equations, dg::DG, cache)
   # to avoid a division by zero if the speed vanishes everywhere,
   # e.g. for steady-state linear advection
   max_scaled_speed = nextfloat(zero(t))
@@ -21,8 +21,8 @@ function max_dt(u::AbstractArray{<:Any,4}, t, mesh::TreeMesh{2},
 end
 
 
-function max_dt(u::AbstractArray{<:Any,4}, t, mesh::TreeMesh{2},
-                constant_speed::Val{true}, equations, dg::DG, cache)
+function max_dt_hyperbolic(u::AbstractArray{<:Any,4}, t, mesh::TreeMesh{2},
+                           constant_speed::Val{true}, equations, dg::DG, cache)
   # to avoid a division by zero if the speed vanishes everywhere,
   # e.g. for steady-state linear advection
   max_scaled_speed = nextfloat(zero(t))
@@ -37,13 +37,13 @@ function max_dt(u::AbstractArray{<:Any,4}, t, mesh::TreeMesh{2},
 end
 
 
-function max_dt(u::AbstractArray{<:Any,4}, t, mesh::ParallelTreeMesh{2},
-                constant_speed::Val{false}, equations, dg::DG, cache)
+function max_dt_hyperbolic(u::AbstractArray{<:Any,4}, t, mesh::ParallelTreeMesh{2},
+                           constant_speed::Val{false}, equations, dg::DG, cache)
   # call the method accepting a general `mesh::TreeMesh{2}`
   # TODO: MPI, we should improve this; maybe we should dispatch on `u`
   #       and create some MPI array type, overloading broadcasting and mapreduce etc.
   #       Then, this specific array type should also work well with DiffEq etc.
-  dt = invoke(max_dt,
+  dt = invoke(max_dt_hyperbolic,
     Tuple{typeof(u), typeof(t), TreeMesh{2},
           typeof(constant_speed), typeof(equations), typeof(dg), typeof(cache)},
     u, t, mesh, constant_speed, equations, dg, cache)
@@ -53,13 +53,13 @@ function max_dt(u::AbstractArray{<:Any,4}, t, mesh::ParallelTreeMesh{2},
 end
 
 
-function max_dt(u::AbstractArray{<:Any,4}, t, mesh::ParallelTreeMesh{2},
-                constant_speed::Val{true}, equations, dg::DG, cache)
+function max_dt_hyperbolic(u::AbstractArray{<:Any,4}, t, mesh::ParallelTreeMesh{2},
+                           constant_speed::Val{true}, equations, dg::DG, cache)
   # call the method accepting a general `mesh::TreeMesh{2}`
   # TODO: MPI, we should improve this; maybe we should dispatch on `u`
   #       and create some MPI array type, overloading broadcasting and mapreduce etc.
   #       Then, this specific array type should also work well with DiffEq etc.
-  dt = invoke(max_dt,
+  dt = invoke(max_dt_hyperbolic,
     Tuple{typeof(u), typeof(t), TreeMesh{2},
           typeof(constant_speed), typeof(equations), typeof(dg), typeof(cache)},
     u, t, mesh, constant_speed, equations, dg, cache)
@@ -68,9 +68,8 @@ function max_dt(u::AbstractArray{<:Any,4}, t, mesh::ParallelTreeMesh{2},
   return dt
 end
 
-
-function max_dt(u::AbstractArray{<:Any,4}, t, mesh::StructuredMesh,
-                constant_speed::Val{false}, equations, dg::DG, cache)
+function max_dt_hyperbolic(u::AbstractArray{<:Any,4}, t, mesh::StructuredMesh,
+                           constant_speed::Val{false}, equations, dg::DG, cache)
   # to avoid a division by zero if the speed vanishes everywhere,
   # e.g. for steady-state linear advection
   max_scaled_speed = nextfloat(zero(t))
@@ -98,8 +97,8 @@ function max_dt(u::AbstractArray{<:Any,4}, t, mesh::StructuredMesh,
 end
 
 
-function max_dt(u::AbstractArray{<:Any,4}, t, mesh::StructuredMesh,
-                constant_speed::Val{true}, equations, dg::DG, cache)
+function max_dt_hyperbolic(u::AbstractArray{<:Any,4}, t, mesh::StructuredMesh,
+                           constant_speed::Val{true}, equations, dg::DG, cache)
   # to avoid a division by zero if the speed vanishes everywhere,
   # e.g. for steady-state linear advection
   max_scaled_speed = nextfloat(zero(t))
@@ -118,3 +117,19 @@ function max_dt(u::AbstractArray{<:Any,4}, t, mesh::StructuredMesh,
 
   return 2 / (nnodes(dg) * max_scaled_speed)
 end
+
+
+function max_dt_parabolic(u::AbstractArray{<:Any,4}, t, mesh::TreeMesh{2},
+                          constant_speed::Val{true}, equations, dg::DG, cache)
+  # to avoid a division by zero if the diffusion vanishes everywhere,
+  # e.g. for steady-state linear advection
+  max_scaled_diffusion = nextfloat(zero(t))
+
+  for element in eachelement(dg, cache)
+    max_λ1, max_λ2 = max_abs_diffusions(equations)
+    inv_jacobian = cache.elements.inverse_jacobian[element]
+    max_scaled_diffusion = max(max_scaled_diffusion, inv_jacobian^2 * (max_λ1 + max_λ2))
+  end
+
+  return 2^2 / (nnodes(dg)^2 * max_scaled_diffusion)
+end
diff --git a/src/callbacks_step/stepsize_dg3d.jl b/src/callbacks_step/stepsize_dg3d.jl
index 38bad06af2e..86661316373 100644
--- a/src/callbacks_step/stepsize_dg3d.jl
+++ b/src/callbacks_step/stepsize_dg3d.jl
@@ -1,6 +1,6 @@
 
-function max_dt(u::AbstractArray{<:Any,5}, t, mesh::TreeMesh{3},
-                constant_speed::Val{false}, equations, dg::DG, cache)
+function max_dt_hyperbolic(u::AbstractArray{<:Any,5}, t, mesh::TreeMesh{3},
+                           constant_speed::Val{false}, equations, dg::DG, cache)
   # to avoid a division by zero if the speed vanishes everywhere,
   # e.g. for steady-state linear advection
   max_scaled_speed = nextfloat(zero(t))
@@ -22,8 +22,8 @@ function max_dt(u::AbstractArray{<:Any,5}, t, mesh::TreeMesh{3},
 end
 
 
-function max_dt(u::AbstractArray{<:Any,5}, t, mesh::TreeMesh{3},
-                constant_speed::Val{true}, equations, dg::DG, cache)
+function max_dt_hyperbolic(u::AbstractArray{<:Any,5}, t, mesh::TreeMesh{3},
+                           constant_speed::Val{true}, equations, dg::DG, cache)
   # to avoid a division by zero if the speed vanishes everywhere,
   # e.g. for steady-state linear advection
   max_scaled_speed = nextfloat(zero(t))
diff --git a/src/equations/equations.jl b/src/equations/equations.jl
index 3b20607e357..1a813f02764 100644
--- a/src/equations/equations.jl
+++ b/src/equations/equations.jl
@@ -192,3 +192,15 @@ include("lattice_boltzmann_3d.jl")
 # Acoustic perturbation equations
 abstract type AbstractAcousticPerturbationEquations{NDIMS, NVARS} <: AbstractEquations{NDIMS, NVARS} end
 include("acoustic_perturbation_2d.jl")
+
+# Gradient equations
+abstract type AbstractGradientEquations{NDIMS, NVARS} <: AbstractEquations{NDIMS, NVARS} end
+include("gradient_equations_2d.jl")
+
+# Heat equation
+abstract type AbstractHeatEquation{NDIMS, NVARS} <: AbstractEquations{NDIMS, NVARS} end
+include("heat_equation_2d.jl")
+
+# Linear scalar advection-diffusion equation
+abstract type AbstractLinearAdvectionDiffusionEquation{NDIMS, NVARS} <: AbstractEquations{NDIMS, NVARS} end
+include("linear_advection_diffusion_2d.jl")
diff --git a/src/equations/gradient_equations_2d.jl b/src/equations/gradient_equations_2d.jl
new file mode 100644
index 00000000000..83639a7769a
--- /dev/null
+++ b/src/equations/gradient_equations_2d.jl
@@ -0,0 +1,65 @@
+
+@doc raw"""
+    GradientEquations2D
+
+The gradient equations
+```math
+q^d - \partial_d u = 0
+```
+in direction `d` in two space dimensions as required for, e.g., the Bassi & Rebay 1 (BR1) or the
+local discontinuous Galerkin (LDG) schemes.
+
+!!! warning "Experimental code"
+    This system of equations is experimental and can change any time.
+"""
+struct GradientEquations2D{RealT<:Real, NVARS} <: AbstractGradientEquations{2, NVARS}
+  orientation::Int
+end
+
+GradientEquations2D(::Type{RealT}, nvars, orientation) where RealT = GradientEquations2D{RealT, nvars}(orientation)
+GradientEquations2D(nvars, orientation) = GradientEquations2D(Float64, nvars, orientation)
+
+
+get_name(::GradientEquations2D) = "GradientEquations2D"
+varnames(::typeof(cons2cons), equations::GradientEquations2D) = SVector(ntuple(v -> "gradient_"*string(v), nvariables(equations)))
+varnames(::typeof(cons2prim), equations::GradientEquations2D) = varnames(cons2cons, equations)
+
+# Set initial conditions at physical location `x` for time `t`
+"""
+    initial_condition_constant(x, t, equations::GradientEquations2D)
+
+A constant initial condition (zero).
+"""
+function initial_condition_constant(x, t, equations::GradientEquations2D)
+  return SVector(ntuple(v -> zero(eltype(x)), nvariables(equations)))
+end
+
+
+function boundary_condition_sin_x(u_inner, orientation, direction, x, t,
+                                  surface_flux_function,
+                                  equations::GradientEquations2D)
+  if orientation == equations.orientation
+    # u_boundary = initial_condition_sin_x(x, t, equation)
+    nu = 1.2e-2
+    omega = pi
+    scalar = sin(omega * x[1]) * exp(-nu * omega^2 * t)
+    u_boundary = SVector(scalar)
+
+    return -u_boundary
+  else
+    return SVector(ntuple(v -> zero(eltype(u_inner)), nvariables(equations)))
+  end
+
+  return flux
+end
+
+# Calculate 1D flux in for a single point
+@inline function calcflux(u, orientation, equations::GradientEquations2D)
+  if orientation == equations.orientation
+    # In the direction specified in the constructor, the flux is equal to the current state.
+    return -u
+  else
+    # Otherwise, the flux is zero.
+    return SVector(ntuple(v -> zero(eltype(u)), nvariables(equations)))
+  end
+end
diff --git a/src/equations/heat_equation_2d.jl b/src/equations/heat_equation_2d.jl
new file mode 100644
index 00000000000..1c15db699f4
--- /dev/null
+++ b/src/equations/heat_equation_2d.jl
@@ -0,0 +1,132 @@
+
+@doc raw"""
+    HeatEquation2D
+
+The heat equation
+```math
+\partial_t u - \nu (\partial^2_1 u + \partial^2_2 u) = 0
+```
+in two space dimensions with constant viscosity ``\nu``.
+
+!!! warning "Experimental code"
+    This system of equations is experimental and can change any time.
+"""
+struct HeatEquation2D{RealT<:Real} <: AbstractHeatEquation{2, 1}
+  nu::RealT
+end
+
+get_name(::HeatEquation2D) = "HeatEquation2D"
+varnames(::typeof(cons2cons), ::HeatEquation2D) = SVector("scalar")
+varnames(::typeof(cons2prim), ::HeatEquation2D) = SVector("scalar")
+
+
+# Set initial conditions at physical location `x` for time `t`
+"""
+    initial_condition_constant(x, t, equations::HeatEquation2D)
+
+A constant initial condition to test free-stream preservation.
+"""
+function initial_condition_constant(x, t, equation::HeatEquation2D)
+  return SVector(2.0)
+end
+
+
+"""
+    initial_condition_convergence_test(x, t, equations::HeatEquation2D)
+
+A smooth initial condition used for convergence tests.
+"""
+function initial_condition_convergence_test(x, t, equation::HeatEquation2D)
+  @unpack nu = equation
+  c = 1.0
+  A = 0.5
+  L = 2
+  f = 1/L
+  omega = 2 * pi * f
+  scalar = c + A * sin(omega * sum(x)) * exp(-2 * nu * omega^2 * t)
+  return SVector(scalar)
+end
+
+
+"""
+    initial_condition_gauss(x, t, equation::HeatEquation2D)
+
+A Gaussian pulse used together with
+[`boundary_condition_gauss`](@ref).
+"""
+function initial_condition_gauss(x, t, equation::HeatEquation2D)
+  return SVector(exp(-(x[1]^2 + x[2]^2)))
+end
+
+
+function initial_condition_sin_x(x, t, equation::HeatEquation2D)
+  @unpack nu = equation
+
+  omega = pi
+  scalar = sin(omega * x[1]) * exp(-nu * omega^2 * t)
+
+  return SVector(scalar)
+end
+
+
+function boundary_condition_sin_x(u_inner, gradients_inner, orientation, direction, x, t,
+                                  surface_flux_function,
+                                  equation::HeatEquation2D)
+  # OBS! This is specifically implemented for BR1
+  return calcflux(u_inner, gradients_inner, orientation, equation)
+end
+
+
+function initial_condition_poisson_periodic(x, t, equation::HeatEquation2D)
+  # elliptic equation: -νΔϕ = f
+  # depending on initial constant state, c, for phi this converges to the solution ϕ + c
+  @unpack nu = equation
+
+  phi = sin(2.0*pi*x[1])*sin(2.0*pi*x[2])*(1 - exp(-8*nu*pi^2*t))
+  return SVector(phi)
+end
+
+
+@inline function source_terms_poisson_periodic(u, x, t, equation::HeatEquation2D)
+  # elliptic equation: -νΔϕ = f
+  # analytical solution: phi = sin(2πx)*sin(2πy) and f = -8νπ^2 sin(2πx)*sin(2πy)
+  C = -8 * equation.nu * pi^2
+
+  x1, x2 = x
+  tmp1 = sinpi(2 * x1)
+  tmp2 = sinpi(2 * x2)
+  du1 = -C*tmp1*tmp2
+
+  return SVector(du1)
+end
+
+
+# Calculate parabolic 1D flux in axis `orientation` for a single point
+@inline function calcflux(u, gradients, orientation, equation::HeatEquation2D)
+  return -equation.nu*gradients[orientation]
+end
+
+
+@inline have_constant_diffusion(::HeatEquation2D) = Val(true)
+
+# FIXME: Find a better name than `max_abs_diffusions` or `max_abs_speeds_viscous`
+@inline function max_abs_diffusions(equation::HeatEquation2D)
+  @unpack nu = equation
+  return (nu, nu)
+end
+
+# Convert conservative variables to primitive
+@inline cons2prim(u, equation::HeatEquation2D) = u
+
+# Convert conservative variables to entropy variables
+@inline cons2entropy(u, equation::HeatEquation2D) = u
+
+
+# Calculate entropy for a conservative state `cons`
+@inline entropy(u::Real, ::HeatEquation2D) = 0.5 * u^2
+@inline entropy(u, equation::HeatEquation2D) = entropy(u[1], equation)
+
+
+# Calculate total energy for a conservative state `cons`
+@inline energy_total(u::Real, ::HeatEquation2D) = 0.5 * u^2
+@inline energy_total(u, equation::HeatEquation2D) = energy_total(u[1], equation)
diff --git a/src/equations/linear_advection_diffusion_2d.jl b/src/equations/linear_advection_diffusion_2d.jl
new file mode 100644
index 00000000000..d2474fab836
--- /dev/null
+++ b/src/equations/linear_advection_diffusion_2d.jl
@@ -0,0 +1,126 @@
+
+@doc raw"""
+    LinearAdvectionDiffusionEquation2D
+
+The linear advection-diffusion equation
+```math
+\partial_t u + a_1 \partial_1 u + a_2 \partial_2 u - \nu (\partial^2_1 u + \partial^2_2 u)= 0
+```
+in two space dimensions with constant velocity `a` and constant viscosity ``\nu``.
+
+!!! warning "Experimental code"
+    This system of equations is experimental and can change any time.
+"""
+struct LinearAdvectionDiffusionEquation2D{RealT<:Real} <: AbstractLinearAdvectionDiffusionEquation{2, 1}
+  advectionvelocity::SVector{2, RealT}
+  nu::RealT
+end
+
+function LinearAdvectionDiffusionEquation2D(a::NTuple{2,<:Real}, nu::Real)
+  LinearAdvectionDiffusionEquation2D(SVector(a), nu)
+end
+
+function LinearAdvectionDiffusionEquation2D(a1::Real, a2::Real, nu::Real)
+  a1, a2, nu = promote(a1, a2, nu)
+  LinearAdvectionDiffusionEquation2D(SVector(a1, a2), nu)
+end
+
+get_name(::LinearAdvectionDiffusionEquation2D) = "LinearAdvectionDiffusionEquation2D"
+varnames(::typeof(cons2cons), ::LinearAdvectionDiffusionEquation2D) = SVector("scalar")
+varnames(::typeof(cons2prim), ::LinearAdvectionDiffusionEquation2D) = SVector("scalar")
+
+
+# Set initial conditions at physical location `x` for time `t`
+"""
+    initial_condition_constant(x, t, equations::LinearAdvectionDiffusionEquation2D)
+
+A constant initial condition to test free-stream preservation.
+"""
+function initial_condition_constant(x, t, equation::LinearAdvectionDiffusionEquation2D)
+  return SVector(2.0)
+end
+
+
+"""
+    initial_condition_convergence_test(x, t, equations::LinearAdvectionDiffusionEquation2D)
+
+A smooth initial condition used for convergence tests.
+"""
+function initial_condition_convergence_test(x, t, equation::LinearAdvectionDiffusionEquation2D)
+  # Store translated coordinate for easy use of exact solution
+  x_trans = x - equation.advectionvelocity * t
+
+  @unpack nu = equation
+  c = 1.0
+  A = 0.5
+  L = 2
+  f = 1/L
+  omega = 2 * pi * f
+  scalar = c + A * sin(omega * sum(x_trans)) * exp(-2 * nu * omega^2 * t)
+  return SVector(scalar)
+end
+
+
+"""
+    initial_condition_gauss(x, t, equation::LinearAdvectionDiffusionEquation2D)
+
+A Gaussian pulse used together with
+[`boundary_condition_gauss`](@ref).
+"""
+function initial_condition_gauss(x, t, equation::LinearAdvectionDiffusionEquation2D)
+  return SVector(exp(-(x[1]^2 + x[2]^2)))
+end
+
+
+# Pre-defined source terms should be implemented as
+# function source_terms_WHATEVER(u, x, t, equations::LinearAdvectionDiffusionEquation2D)
+
+
+# Calculate hyperbolic 1D flux in axis `orientation` for a single point
+@inline function calcflux(u, orientation, equation::LinearAdvectionDiffusionEquation2D)
+  a = equation.advectionvelocity[orientation]
+  return a * u
+end
+
+# Calculate parabolic 1D flux in axis `orientation` for a single point
+@inline function calcflux(u, gradients, orientation, equation::LinearAdvectionDiffusionEquation2D)
+  return -equation.nu*gradients[orientation]
+end
+
+
+function flux_lax_friedrichs(u_ll, u_rr, orientation, equation::LinearAdvectionDiffusionEquation2D)
+  a = equation.advectionvelocity[orientation]
+  return 0.5 * ( a * (u_ll + u_rr) - abs(a) * (u_rr - u_ll) )
+end
+
+
+
+@inline have_constant_speed(::LinearAdvectionDiffusionEquation2D) = Val(true)
+
+@inline function max_abs_speeds(equation::LinearAdvectionDiffusionEquation2D)
+  return abs.(equation.advectionvelocity)
+end
+
+@inline have_constant_diffusion(::LinearAdvectionDiffusionEquation2D) = Val(true)
+
+@inline function max_abs_diffusions(equation::LinearAdvectionDiffusionEquation2D)
+  @unpack nu = equation
+  return (nu, nu)
+end
+
+
+# Convert conservative variables to primitive
+@inline cons2prim(u, equation::LinearAdvectionDiffusionEquation2D) = u
+
+# Convert conservative variables to entropy variables
+@inline cons2entropy(u, equation::LinearAdvectionDiffusionEquation2D) = u
+
+
+# Calculate entropy for a conservative state `cons`
+@inline entropy(u::Real, ::LinearAdvectionDiffusionEquation2D) = 0.5 * u^2
+@inline entropy(u, equation::LinearAdvectionDiffusionEquation2D) = entropy(u[1], equation)
+
+
+# Calculate total energy for a conservative state `cons`
+@inline energy_total(u::Real, ::LinearAdvectionDiffusionEquation2D) = 0.5 * u^2
+@inline energy_total(u, equation::LinearAdvectionDiffusionEquation2D) = energy_total(u[1], equation)
diff --git a/src/semidiscretization/semidiscretization_euler_gravity.jl b/src/semidiscretization/semidiscretization_euler_gravity.jl
index fa5ca70b2b4..a65cc0453ac 100644
--- a/src/semidiscretization/semidiscretization_euler_gravity.jl
+++ b/src/semidiscretization/semidiscretization_euler_gravity.jl
@@ -177,6 +177,17 @@ end
 end
 
 
+function max_dt(u_ode::AbstractVector, t, cfl_number::Real, semi::SemidiscretizationEulerGravity)
+  mesh, equations, solver, cache = mesh_equations_solver_cache(semi)
+  u = wrap_array(u_ode, mesh, equations, solver, cache)
+
+  dt = cfl_number * max_dt_hyperbolic(u, t, mesh, have_constant_speed(equations), equations, solver,
+                                      cache)
+
+  return dt
+end
+
+
 function rhs!(du_ode, u_ode, semi::SemidiscretizationEulerGravity, t)
   @unpack semi_euler, semi_gravity, cache = semi
 
@@ -241,9 +252,12 @@ function update_gravity!(semi::SemidiscretizationEulerGravity, u_ode::AbstractVe
   # iterate gravity solver until convergence or maximum number of iterations are reached
   @unpack equations = semi_gravity
   while !finalstep
-    dt = @timeit_debug timer() "calculate dt" cfl * max_dt(u_gravity, t, semi_gravity.mesh,
-                                                           have_constant_speed(equations), equations,
-                                                           semi_gravity.solver, semi_gravity.cache)
+    dt = @timeit_debug timer() "calculate dt" cfl * max_dt_hyperbolic(u_gravity, t,
+                                                                      semi_gravity.mesh,
+                                                                      have_constant_speed(equations),
+                                                                      equations,
+                                                                      semi_gravity.solver,
+                                                                      semi_gravity.cache)
 
     # evolve solution by one pseudo-time step
     time_start = time_ns()
diff --git a/src/semidiscretization/semidiscretization_hyperbolic.jl b/src/semidiscretization/semidiscretization_hyperbolic.jl
index 39e25454eff..76f4750f8a6 100644
--- a/src/semidiscretization/semidiscretization_hyperbolic.jl
+++ b/src/semidiscretization/semidiscretization_hyperbolic.jl
@@ -191,6 +191,17 @@ function compute_coefficients!(u_ode::AbstractVector, t, semi::Semidiscretizatio
 end
 
 
+function max_dt(u_ode::AbstractVector, t, cfl_number::Real, semi::SemidiscretizationHyperbolic)
+  mesh, equations, solver, cache = mesh_equations_solver_cache(semi)
+  u = wrap_array(u_ode, mesh, equations, solver, cache)
+
+  dt = cfl_number * max_dt_hyperbolic(u, t, mesh, have_constant_speed(equations), equations, solver,
+                                      cache)
+
+  return dt
+end
+
+
 function rhs!(du_ode, u_ode, semi::SemidiscretizationHyperbolic, t)
   @unpack mesh, equations, initial_condition, boundary_conditions, source_terms, solver, cache = semi
 
diff --git a/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl b/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl
new file mode 100644
index 00000000000..a92876a829e
--- /dev/null
+++ b/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl
@@ -0,0 +1,203 @@
+
+# TODO: Taal refactor, Mesh<:AbstractMesh{NDIMS}, Equations<:AbstractEquations{NDIMS} ?
+"""
+    SemidiscretizationHyperbolicParabolic
+
+A struct containing everything needed to describe a spatial semidiscretization
+of a hyperbolic-parabolic conservation law.
+
+!!! warning "Experimental code"
+    This semidiscretization is experimental and can change any time.
+"""
+struct SemidiscretizationHyperbolicParabolic{SemiHyperbolic, SemiParabolic, Cache} <: AbstractSemidiscretization
+  semi_hyperbolic::SemiHyperbolic
+  semi_parabolic::SemiParabolic
+  cache::Cache
+  performance_counter::PerformanceCounter
+
+  function SemidiscretizationHyperbolicParabolic{SemiHyperbolic, SemiParabolic, Cache}(
+      semi_hyperbolic::SemiHyperbolic, semi_parabolic::SemiParabolic,
+      cache::Cache) where {SemiHyperbolic, SemiParabolic, Cache}
+
+    performance_counter = PerformanceCounter()
+
+    new(semi_hyperbolic, semi_parabolic, cache, performance_counter)
+  end
+end
+
+"""
+    SemidiscretizationHyperbolicParabolic(semi_hyperbolic, semi_parabolic)
+
+Construct a semidiscretization of a hyperbolic-parabolic PDE by combining the purely hyperbolic and
+purely parabolic semidiscretizations.
+
+!!! warning "Experimental code"
+    This semidiscretization is experimental and can change any time.
+"""
+function SemidiscretizationHyperbolicParabolic(semi_hyperbolic, semi_parabolic)
+  cache = (; du_ode_parabolic=Vector{real(semi_parabolic)}())
+
+  SemidiscretizationHyperbolicParabolic{typeof(semi_hyperbolic), typeof(semi_parabolic), typeof(cache)}(
+    semi_hyperbolic, semi_parabolic, cache)
+end
+
+
+function Base.show(io::IO, semi::SemidiscretizationHyperbolicParabolic)
+  print(io, "SemidiscretizationHyperbolicParabolic(")
+  print(io,       semi.semi_hyperbolic)
+  print(io, ", ", semi.semi_parabolic)
+  print(io, ", cache(")
+  for (idx,key) in enumerate(keys(semi.cache))
+    idx > 1 && print(io, " ")
+    print(io, key)
+  end
+  print(io, "))")
+end
+
+function Base.show(io::IO, ::MIME"text/plain", semi::SemidiscretizationHyperbolicParabolic)
+  if get(io, :compact, false)
+    show(io, semi)
+  else
+    summary_header(io, "SemidiscretizationHyperbolicParabolic")
+    summary_line(io, "hyperbolic system", semi.semi_hyperbolic)
+    summary_line(io, "parabolic system", semi.semi_parabolic)
+    summary_footer(io)
+  end
+end
+
+
+@inline Base.ndims(semi::SemidiscretizationHyperbolicParabolic) = ndims(semi.semi_hyperbolic.mesh)
+
+@inline nvariables(semi::SemidiscretizationHyperbolicParabolic) = nvariables(semi.semi_hyperbolic.equations)
+
+@inline Base.real(semi::SemidiscretizationHyperbolicParabolic) = real(semi.semi_hyperbolic.solver)
+
+
+@inline function mesh_equations_solver_cache(semi::SemidiscretizationHyperbolicParabolic)
+  return mesh_equations_solver_cache(semi.semi_hyperbolic)
+end
+
+
+function calc_error_norms(func, u_ode::AbstractVector, t, analyzer, semi::SemidiscretizationHyperbolicParabolic, cache_analysis)
+  mesh, equations, solver, cache = mesh_equations_solver_cache(semi.semi_hyperbolic)
+  @unpack initial_condition = semi.semi_hyperbolic
+
+  u = wrap_array(u_ode, mesh, equations, solver, cache)
+
+  calc_error_norms(func, u, t, analyzer, mesh, equations, initial_condition, solver, cache, cache_analysis)
+end
+
+function calc_error_norms(func, u, t, analyzer, semi::SemidiscretizationHyperbolicParabolic, cache_analysis)
+  mesh, equations, solver, cache = mesh_equations_solver_cache(semi.semi_hyperbolic)
+  @unpack initial_condition = semi.semi_hyperbolic
+
+  calc_error_norms(func, u, t, analyzer, mesh, equations, initial_condition, solver, cache, cache_analysis)
+end
+
+
+function compute_coefficients(t, semi::SemidiscretizationHyperbolicParabolic)
+  compute_coefficients(semi.semi_hyperbolic.initial_condition, t, semi.semi_hyperbolic)
+end
+
+function compute_coefficients!(u_ode::AbstractVector, t, semi::SemidiscretizationHyperbolicParabolic)
+  compute_coefficients!(u_ode, semi.semi_hyperbolic.initial_condition, t, semi.semi_hyperbolic)
+end
+
+
+function max_dt(u_ode::AbstractVector, t, cfl_number,
+                semi::SemidiscretizationHyperbolicParabolic)
+  @unpack semi_hyperbolic, semi_parabolic = semi
+  dt = min(max_dt(u_ode, t, cfl_number[1], semi_hyperbolic),
+           max_dt(u_ode, t, cfl_number[2], semi_parabolic))
+
+  return dt
+end
+
+
+function rhs!(du_ode, u_ode, semi::SemidiscretizationHyperbolicParabolic, t)
+  @unpack du_ode_parabolic = semi.cache
+
+  # TODO: Taal decide, do we need to pass the mesh?
+  time_start = time_ns()
+
+  # Resize the storage for the parabolic time derivative in case it was changed
+  resize!(du_ode_parabolic, length(du_ode))
+
+  @timeit_debug timer() "rhs!" begin
+    # Compute hyperbolic time derivative
+    rhs!(du_ode, u_ode, semi.semi_hyperbolic, t)
+
+    # Compute parabolic time derivative
+    rhs!(du_ode_parabolic, u_ode, semi.semi_parabolic, t)
+
+    # Add parabolic update to hyperbolic time derivative
+    du_ode .+= du_ode_parabolic
+  end
+
+  runtime = time_ns() - time_start
+  put!(semi.performance_counter, runtime)
+
+  return nothing
+end
+
+
+"""
+    SemidiscretizationHyperbolicParabolicBR1(mesh, equations, initial_condition, solver::DGSEM;
+                                             source_terms=nothing,
+                                             boundary_conditions=boundary_condition_periodic)
+
+Convenience function to construct a semidiscretization of a hyperbolic-parabolic PDE with the BR1
+scheme for a DGSEM solver. The passed arguments must correspond to the hyperbolic setup, including
+the `solver`. Internally, a matching solver for the parabolic system is created and both a
+hyperbolic and a parabolic semidiscretization are created and passed to a
+`SemidiscretizationHyperbolicParabolic`.
+
+!!! warning "Experimental code"
+    This semidiscretization is experimental and can change any time.
+"""
+function SemidiscretizationHyperbolicParabolicBR1(mesh, equations, initial_condition, solver::DGSEM;
+                                                  source_terms=nothing,
+                                                  boundary_conditions=boundary_condition_periodic,
+                                                  RealT=real(solver),
+                                                  volume_integral_parabolic=VolumeIntegralWeakForm())
+  semi_hyperbolic = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver;
+                                                 source_terms=source_terms,
+                                                 boundary_conditions=boundary_conditions,
+                                                 RealT=RealT)
+
+  solver_parabolic = DGSEM(solver.basis, flux_central, volume_integral_parabolic, solver.mortar)
+  semi_parabolic = SemidiscretizationParabolicAuxVars(mesh, equations, initial_condition,
+                                                      solver_parabolic,
+                                                      source_terms=source_terms,
+                                                      boundary_conditions=boundary_conditions,
+                                                      RealT=RealT)
+
+  return SemidiscretizationHyperbolicParabolic(semi_hyperbolic, semi_parabolic)
+end
+
+
+@inline function (amr_callback::AMRCallback)(u_ode::AbstractVector,
+                                             semi::SemidiscretizationHyperbolicParabolic,
+                                             t, iter; kwargs...)
+  # Get semidiscretizations
+  @unpack semi_hyperbolic, semi_parabolic = semi
+  @unpack semi_gradients_x, semi_gradients_y = semi_parabolic.cache
+
+  # Perform AMR based on hyperbolic solver
+  has_changed = amr_callback(u_ode, mesh_equations_solver_cache(semi_hyperbolic)..., t, iter;
+                             kwargs...)
+
+  if has_changed
+    # Reinitialize data structures for all secondary solvers (no solution data conversion needed)
+    reinitialize_containers!(mesh_equations_solver_cache(semi_parabolic)...)
+    reinitialize_containers!(mesh_equations_solver_cache(semi_gradients_x)...)
+    reinitialize_containers!(mesh_equations_solver_cache(semi_gradients_y)...)
+
+    # Resize storage for auxiliary variables (= gradients)
+    @unpack u_ode_gradients_x, u_ode_gradients_y = semi_parabolic.cache
+    resize!(u_ode_gradients_x, length(u_ode))
+    resize!(u_ode_gradients_y, length(u_ode))
+  end
+
+  return has_changed
+end
diff --git a/src/semidiscretization/semidiscretization_parabolic_auxvars.jl b/src/semidiscretization/semidiscretization_parabolic_auxvars.jl
new file mode 100644
index 00000000000..0d7f8a53543
--- /dev/null
+++ b/src/semidiscretization/semidiscretization_parabolic_auxvars.jl
@@ -0,0 +1,224 @@
+
+# TODO: Taal refactor, Mesh<:AbstractMesh{NDIMS}, Equations<:AbstractEquations{NDIMS} ?
+"""
+    SemidiscretizationParabolicAuxVars
+
+A struct containing everything needed to describe a spatial semidiscretization
+of a parabolic conservation law.
+
+!!! warning "Experimental code"
+    This semidiscretization is experimental and can change any time.
+"""
+struct SemidiscretizationParabolicAuxVars{Mesh, Equations, InitialCondition, BoundaryConditions,
+                                          SourceTerms, Solver, Cache} <: AbstractSemidiscretization
+
+  mesh::Mesh
+  equations::Equations
+
+  # This guy is a bit messy since we abuse it as some kind of "exact solution"
+  # although this doesn't really exist...
+  initial_condition::InitialCondition
+
+  boundary_conditions::BoundaryConditions
+  source_terms::SourceTerms
+  solver::Solver
+  cache::Cache
+  performance_counter::PerformanceCounter
+
+  function SemidiscretizationParabolicAuxVars{Mesh, Equations, InitialCondition, BoundaryConditions,
+                                              SourceTerms, Solver, Cache}(
+      mesh::Mesh, equations::Equations,
+      initial_condition::InitialCondition, boundary_conditions::BoundaryConditions,
+      source_terms::SourceTerms,
+      solver::Solver, cache::Cache) where {Mesh, Equations, InitialCondition, BoundaryConditions, SourceTerms, Solver, Cache}
+    @assert ndims(mesh) == ndims(equations)
+
+    performance_counter = PerformanceCounter()
+
+    new(mesh, equations, initial_condition, boundary_conditions, source_terms, solver, cache, performance_counter)
+  end
+end
+
+"""
+    SemidiscretizationParabolicAuxVars(mesh, equations, initial_condition, solver;
+                                       source_terms=nothing,
+                                       boundary_conditions=boundary_condition_periodic)
+
+Construct a semidiscretization of a parabolic PDE.
+
+!!! warning "Experimental code"
+    This semidiscretization is experimental and can change any time.
+"""
+function SemidiscretizationParabolicAuxVars(mesh, equations, initial_condition, solver;
+                                            source_terms=nothing,
+                                            boundary_conditions=boundary_condition_periodic,
+                                            RealT=real(solver))
+
+  cache = create_cache(mesh, equations, solver, RealT)
+  _boundary_conditions = digest_boundary_conditions(boundary_conditions)
+
+  solver_auxvars = create_solver_auxvars(solver)
+
+  equations_gradients_x = GradientEquations2D(nvariables(equations), 1)
+  semi_gradients_x = SemidiscretizationHyperbolic(mesh, equations_gradients_x,
+                                                  initial_condition_constant, solver_auxvars,
+                                                  boundary_conditions=boundary_conditions)
+
+  equations_gradients_y = GradientEquations2D(nvariables(equations), 2)
+  semi_gradients_y = SemidiscretizationHyperbolic(mesh, equations_gradients_y,
+                                                  initial_condition_constant, solver_auxvars,
+                                                  boundary_conditions=boundary_conditions)
+
+  u_ode_gradients_x = compute_coefficients(nothing, semi_gradients_x)
+  u_ode_gradients_y = compute_coefficients(nothing, semi_gradients_y)
+
+  cache = (; cache..., semi_gradients_x, u_ode_gradients_x, semi_gradients_y, u_ode_gradients_y)
+
+  SemidiscretizationParabolicAuxVars{typeof(mesh), typeof(equations), typeof(initial_condition), typeof(_boundary_conditions), typeof(source_terms), typeof(solver), typeof(cache)}(
+    mesh, equations, initial_condition, _boundary_conditions, source_terms, solver, cache)
+end
+
+
+function Base.show(io::IO, semi::SemidiscretizationParabolicAuxVars)
+  print(io, "SemidiscretizationParabolicAuxVars(")
+  print(io,       semi.mesh)
+  print(io, ", ", semi.equations)
+  print(io, ", ", semi.initial_condition)
+  print(io, ", ", semi.boundary_conditions)
+  print(io, ", ", semi.source_terms)
+  print(io, ", ", semi.solver)
+  print(io, ", cache(")
+  for (idx,key) in enumerate(keys(semi.cache))
+    idx > 1 && print(io, " ")
+    print(io, key)
+  end
+  print(io, "))")
+end
+
+function Base.show(io::IO, ::MIME"text/plain", semi::SemidiscretizationParabolicAuxVars)
+  if get(io, :compact, false)
+    show(io, semi)
+  else
+    summary_header(io, "SemidiscretizationParabolicAuxVars")
+    summary_line(io, "#spatial dimensions", ndims(semi.equations))
+    summary_line(io, "mesh", semi.mesh)
+    summary_line(io, "equations", typeof(semi.equations).name)
+    summary_line(io, "initial condition", semi.initial_condition)
+    summary_line(io, "boundary conditions", 2*ndims(semi))
+    if (semi.boundary_conditions isa Tuple ||
+        semi.boundary_conditions isa NamedTuple ||
+        semi.boundary_conditions isa AbstractArray)
+      bcs = semi.boundary_conditions
+    else
+      bcs = collect(semi.boundary_conditions for _ in 1:(2*ndims(semi)))
+    end
+    summary_line(increment_indent(io), "negative x", bcs[1])
+    summary_line(increment_indent(io), "positive x", bcs[2])
+    if ndims(semi) > 1
+      summary_line(increment_indent(io), "negative y", bcs[3])
+      summary_line(increment_indent(io), "positive y", bcs[4])
+    end
+    if ndims(semi) > 2
+      summary_line(increment_indent(io), "negative z", bcs[5])
+      summary_line(increment_indent(io), "positive z", bcs[6])
+    end
+    summary_line(io, "source terms", semi.source_terms)
+    summary_line(io, "solver", typeof(semi.solver).name)
+    summary_line(io, "total #DOFs", ndofs(semi))
+    summary_footer(io)
+  end
+end
+
+
+@inline Base.ndims(semi::SemidiscretizationParabolicAuxVars) = ndims(semi.mesh)
+
+@inline nvariables(semi::SemidiscretizationParabolicAuxVars) = nvariables(semi.equations)
+
+@inline Base.real(semi::SemidiscretizationParabolicAuxVars) = real(semi.solver)
+
+
+@inline function mesh_equations_solver_cache(semi::SemidiscretizationParabolicAuxVars)
+  @unpack mesh, equations, solver, cache = semi
+  return mesh, equations, solver, cache
+end
+
+
+function calc_error_norms(func, u_ode::AbstractVector, t, analyzer, semi::SemidiscretizationParabolicAuxVars, cache_analysis)
+  @unpack mesh, equations, initial_condition, solver, cache = semi
+  u = wrap_array(u_ode, mesh, equations, solver, cache)
+
+  calc_error_norms(func, u, t, analyzer, mesh, equations, initial_condition, solver, cache, cache_analysis)
+end
+
+function calc_error_norms(func, u, t, analyzer, semi::SemidiscretizationParabolicAuxVars, cache_analysis)
+  @unpack mesh, equations, initial_condition, solver, cache = semi
+
+  calc_error_norms(func, u, t, analyzer, mesh, equations, initial_condition, solver, cache, cache_analysis)
+end
+
+
+function compute_coefficients(t, semi::SemidiscretizationParabolicAuxVars)
+  compute_coefficients(semi.initial_condition, t, semi)
+end
+
+function compute_coefficients!(u_ode::AbstractVector, t, semi::SemidiscretizationParabolicAuxVars)
+  compute_coefficients!(u_ode, semi.initial_condition, t, semi)
+end
+
+
+function max_dt(u_ode::AbstractVector, t, cfl_number::Real, semi::SemidiscretizationParabolicAuxVars)
+  mesh, equations, solver, cache = mesh_equations_solver_cache(semi)
+  u = wrap_array(u_ode, mesh, equations, solver, cache)
+
+  dt = cfl_number * max_dt_parabolic(u, t, mesh, have_constant_diffusion(equations), equations,
+                                     solver, cache)
+
+  return dt
+end
+
+
+function rhs!(du_ode, u_ode, semi::SemidiscretizationParabolicAuxVars, t)
+  @unpack mesh, equations, initial_condition, boundary_conditions, source_terms, solver, cache = semi
+
+  u  = wrap_array(u_ode,  mesh, equations, solver, cache)
+  du = wrap_array(du_ode, mesh, equations, solver, cache)
+
+  # Compute gradients
+  @unpack semi_gradients_x, u_ode_gradients_x, semi_gradients_y, u_ode_gradients_y = cache
+  rhs!(u_ode_gradients_x, u_ode, semi_gradients_x, t)
+  rhs!(u_ode_gradients_y, u_ode, semi_gradients_y, t)
+  gradients = (wrap_array(u_ode_gradients_x, mesh, equations, solver, cache),
+               wrap_array(u_ode_gradients_y, mesh, equations, solver, cache))
+  cache_gradients = (semi_gradients_x.cache, semi_gradients_y.cache)
+
+  # TODO: Taal decide, do we need to pass the mesh?
+  time_start = time_ns()
+  @timeit_debug timer() "rhs!" rhs!(du, u, gradients, t, mesh, equations, initial_condition, boundary_conditions, source_terms, solver, cache, cache_gradients)
+  runtime = time_ns() - time_start
+  put!(semi.performance_counter, runtime)
+
+  return nothing
+end
+
+
+@inline function (amr_callback::AMRCallback)(u_ode::AbstractVector,
+                                             semi::SemidiscretizationParabolicAuxVars,
+                                             t, iter;
+                                             kwargs...)
+  # Perform AMR based on main solver
+  has_changed = amr_callback(u_ode, mesh_equations_solver_cache(semi)..., t, iter; kwargs...)
+
+  if has_changed
+    # Reinitialize data structures for all secondary solvers (no solution data conversion needed)
+    @unpack semi_gradients_x, semi_gradients_y = semi.cache
+    reinitialize_containers!(mesh_equations_solver_cache(semi_gradients_x)...)
+    reinitialize_containers!(mesh_equations_solver_cache(semi_gradients_y)...)
+
+    # Resize storage for auxiliary variables (= gradients)
+    @unpack u_ode_gradients_x, u_ode_gradients_y = semi.cache
+    resize!(u_ode_gradients_x, length(u_ode))
+    resize!(u_ode_gradients_y, length(u_ode))
+  end
+
+  return has_changed
+end
diff --git a/src/solvers/dg/dg.jl b/src/solvers/dg/dg.jl
index 27356343322..5251ae64e27 100644
--- a/src/solvers/dg/dg.jl
+++ b/src/solvers/dg/dg.jl
@@ -321,6 +321,13 @@ function volume_jacobian(element, mesh::TreeMesh, cache)
 end
 
 
+function create_solver_auxvars(dg::DGSEM)
+  @unpack basis, mortar = dg
+
+  return DGSEM(basis, flux_central, VolumeIntegralWeakForm(), mortar)
+end
+
+
 
 # indicators used for shock-capturing and AMR
 include("indicators.jl")
@@ -336,6 +343,7 @@ include("dg_1d.jl")
 include("containers_2d.jl")
 include("dg_2d.jl")
 include("dg_2d_parallel.jl")
+include("dg_2d_parabolic.jl")
 
 # 3D DG implementation
 include("containers_3d.jl")
diff --git a/src/solvers/dg/dg_2d_parabolic.jl b/src/solvers/dg/dg_2d_parabolic.jl
new file mode 100644
index 00000000000..9dbf0b6f5d4
--- /dev/null
+++ b/src/solvers/dg/dg_2d_parabolic.jl
@@ -0,0 +1,287 @@
+# This file collects all methods that have been updated to work with parabolic systems of equations
+
+function rhs!(du::AbstractArray{<:Any,4}, u, gradients, t,
+              mesh::TreeMesh{2}, equations,
+              initial_condition, boundary_conditions, source_terms,
+              dg::DG, cache, cache_gradients)
+  # Reset du
+  @timeit_debug timer() "reset ∂u/∂t" du .= zero(eltype(du))
+
+  # Calculate volume integral
+  @timeit_debug timer() "volume integral" calc_volume_integral!(du, u, gradients, have_nonconservative_terms(equations), equations,
+                                                                dg.volume_integral, dg, cache)
+
+  # Prolong solution and gradients to interfaces
+  @timeit_debug timer() "prolong2interfaces" prolong2interfaces!(cache, cache_gradients, u, gradients, equations, dg)
+
+  # Calculate interface fluxes
+  @timeit_debug timer() "interface flux" calc_interface_flux!(cache.elements.surface_flux_values,
+                                                              have_nonconservative_terms(equations), equations,
+                                                              dg, cache, cache_gradients)
+
+  # Prolong solution and gradients to boundaries
+  @timeit_debug timer() "prolong2boundaries" prolong2boundaries!(cache, cache_gradients, u, gradients, equations, dg)
+
+  # Calculate boundary fluxes
+  @timeit_debug timer() "boundary flux" calc_boundary_flux!(cache, cache_gradients, t, boundary_conditions, equations, dg)
+
+  # Prolong solution and gradients to mortars
+  @timeit_debug timer() "prolong2mortars" prolong2mortars!(cache, cache_gradients, u, gradients,
+                                                           equations, dg.mortar, dg)
+
+  # Calculate mortar fluxes
+  @timeit_debug timer() "mortar flux" calc_mortar_flux!(cache.elements.surface_flux_values,
+                                                        have_nonconservative_terms(equations), equations,
+                                                        dg.mortar, dg, cache, cache_gradients)
+
+  # Calculate surface integrals
+  @timeit_debug timer() "surface integral" calc_surface_integral!(du, equations, dg, cache)
+
+  # Apply Jacobian from mapping to reference element
+  @timeit_debug timer() "Jacobian" apply_jacobian!(du, equations, dg, cache)
+
+  # Calculate source terms
+  @timeit_debug timer() "source terms" calc_sources!(du, u, t, source_terms, equations, dg, cache)
+
+  return nothing
+end
+
+
+function calc_volume_integral!(du::AbstractArray{<:Any,4}, u, gradients,
+                               nonconservative_terms::Val{false}, equations,
+                               volume_integral::VolumeIntegralWeakForm,
+                               dg::DGSEM, cache)
+  @unpack derivative_dhat = dg.basis
+
+  Threads.@threads for element in eachelement(dg, cache)
+    for j in eachnode(dg), i in eachnode(dg)
+      u_node = get_node_vars(u, equations, dg, i, j, element)
+      gradients_x_node = get_node_vars(gradients[1], equations, dg, i, j, element)
+      gradients_y_node = get_node_vars(gradients[2], equations, dg, i, j, element)
+      gradients_node = (gradients_x_node, gradients_y_node)
+
+      flux1 = calcflux(u_node, gradients_node, 1, equations)
+      for ii in eachnode(dg)
+        integral_contribution = derivative_dhat[ii, i] * flux1
+        add_to_node_vars!(du, integral_contribution, equations, dg, ii, j, element)
+      end
+
+      flux2 = calcflux(u_node, gradients_node, 2, equations)
+      for jj in eachnode(dg)
+        integral_contribution = derivative_dhat[jj, j] * flux2
+        add_to_node_vars!(du, integral_contribution, equations, dg, i, jj, element)
+      end
+    end
+  end
+
+  return nothing
+end
+
+
+function prolong2interfaces!(cache, cache_gradients, u::AbstractArray{<:Any,4}, gradients,
+                             equations, dg::DG)
+  prolong2interfaces!(cache, u, equations, dg)
+  prolong2interfaces!(cache_gradients[1], gradients[1], equations, dg)
+  prolong2interfaces!(cache_gradients[2], gradients[2], equations, dg)
+
+  return nothing
+end
+
+
+function calc_interface_flux!(surface_flux_values::AbstractArray{<:Any,4},
+                              nonconservative_terms::Val{false}, equations,
+                              dg::DG, cache, cache_gradients)
+  @unpack surface_flux = dg
+  @unpack u, neighbor_ids, orientations = cache.interfaces
+  gradients_x = cache_gradients[1].interfaces.u
+  gradients_y = cache_gradients[2].interfaces.u
+
+  Threads.@threads for interface in eachinterface(dg, cache)
+    # Get neighboring elements
+    left_id  = neighbor_ids[1, interface]
+    right_id = neighbor_ids[2, interface]
+
+    # Determine interface direction with respect to elements:
+    # orientation = 1: left -> 2, right -> 1
+    # orientation = 2: left -> 4, right -> 3
+    left_direction  = 2 * orientations[interface]
+    right_direction = 2 * orientations[interface] - 1
+
+    for i in eachnode(dg)
+      # Call pointwise Riemann solver
+      u_ll, u_rr = get_surface_node_vars(u, equations, dg, i, interface)
+      gradients_x_ll, gradients_x_rr = get_surface_node_vars(gradients_x, equations, dg, i, interface)
+      gradients_y_ll, gradients_y_rr = get_surface_node_vars(gradients_y, equations, dg, i, interface)
+      gradients_ll = (gradients_x_ll, gradients_y_ll)
+      gradients_rr = (gradients_x_rr, gradients_y_rr)
+      flux = surface_flux(u_ll, u_rr, gradients_ll, gradients_rr, orientations[interface], equations)
+
+      # Copy flux to left and right element storage
+      for v in eachvariable(equations)
+        surface_flux_values[v, i, left_direction,  left_id]  = flux[v]
+        surface_flux_values[v, i, right_direction, right_id] = flux[v]
+      end
+    end
+  end
+
+  return nothing
+end
+
+
+function prolong2boundaries!(cache, cache_gradients, u::AbstractArray{<:Any,4}, gradients,
+                             equations, dg::DG)
+  prolong2boundaries!(cache, u, equations, dg)
+  prolong2boundaries!(cache_gradients[1], gradients[1], equations, dg)
+  prolong2boundaries!(cache_gradients[2], gradients[2], equations, dg)
+
+  return nothing
+end
+
+
+# TODO: Taal dimension agnostic
+function calc_boundary_flux!(cache, cache_gradients, t,
+                             boundary_condition::BoundaryConditionPeriodic,
+                             equations::AbstractEquations{2}, dg::DG)
+  @assert isempty(eachboundary(dg, cache))
+end
+
+# TODO: Taal dimension agnostic
+function calc_boundary_flux!(cache, cache_gradients, t, boundary_condition,
+                             equations::AbstractEquations{2}, dg::DG)
+  @unpack surface_flux_values = cache.elements
+  @unpack n_boundaries_per_direction = cache.boundaries
+
+  # Calculate indices
+  lasts = accumulate(+, n_boundaries_per_direction)
+  firsts = lasts - n_boundaries_per_direction .+ 1
+
+  # Calc boundary fluxes in each direction
+  for direction in eachindex(firsts)
+    calc_boundary_flux_by_direction!(surface_flux_values, t, boundary_condition,
+                                     equations, dg, cache, cache_gradients,
+                                     direction, firsts[direction], lasts[direction])
+  end
+end
+
+function calc_boundary_flux!(cache, cache_gradients, t,
+                             boundary_conditions::Union{NamedTuple,Tuple},
+                             equations::AbstractEquations{2}, dg::DG)
+  @unpack surface_flux_values = cache.elements
+  @unpack n_boundaries_per_direction = cache.boundaries
+
+  # Calculate indices
+  lasts = accumulate(+, n_boundaries_per_direction)
+  firsts = lasts - n_boundaries_per_direction .+ 1
+
+  # Calc boundary fluxes in each direction
+  calc_boundary_flux_by_direction!(surface_flux_values, t, boundary_conditions[1],
+                                   equations, dg, cache, cache_gradients, 1, firsts[1], lasts[1])
+  calc_boundary_flux_by_direction!(surface_flux_values, t, boundary_conditions[2],
+                                   equations, dg, cache, cache_gradients, 2, firsts[2], lasts[2])
+  calc_boundary_flux_by_direction!(surface_flux_values, t, boundary_conditions[3],
+                                   equations, dg, cache, cache_gradients, 3, firsts[3], lasts[3])
+  calc_boundary_flux_by_direction!(surface_flux_values, t, boundary_conditions[4],
+                                   equations, dg, cache, cache_gradients, 4, firsts[4], lasts[4])
+end
+
+function calc_boundary_flux_by_direction!(surface_flux_values::AbstractArray{<:Any,4}, t,
+                                          boundary_condition, equations, dg::DG,
+                                          cache, cache_gradients,
+                                          direction, first_boundary, last_boundary)
+  @unpack surface_flux = dg
+  @unpack u, neighbor_ids, neighbor_sides, node_coordinates, orientations = cache.boundaries
+  gradients_x = cache_gradients[1].boundaries.u
+  gradients_y = cache_gradients[2].boundaries.u
+
+  Threads.@threads for boundary in first_boundary:last_boundary
+    # Get neighboring element
+    neighbor = neighbor_ids[boundary]
+
+    for i in eachnode(dg)
+      # Get boundary flux
+      u_ll, u_rr = get_surface_node_vars(u, equations, dg, i, boundary)
+      gradients_x_ll, gradients_x_rr = get_surface_node_vars(gradients_x, equations, dg, i, boundary)
+      gradients_y_ll, gradients_y_rr = get_surface_node_vars(gradients_y, equations, dg, i, boundary)
+      gradients_ll = (gradients_x_ll, gradients_y_ll)
+      gradients_rr = (gradients_x_rr, gradients_y_rr)
+      if neighbor_sides[boundary] == 1 # Element is on the left, boundary on the right
+        u_inner = u_ll
+        gradients_inner = gradients_ll
+      else # Element is on the right, boundary on the left
+        u_inner = u_rr
+        gradients_inner = gradients_rr
+      end
+      x = get_node_coords(node_coordinates, equations, dg, i, boundary)
+      flux = boundary_condition(u_inner, gradients_inner, orientations[boundary], direction, x, t,
+                                surface_flux, equations)
+
+      # Copy flux to left and right element storage
+      for v in eachvariable(equations)
+        surface_flux_values[v, i, direction, neighbor] = flux[v]
+      end
+    end
+  end
+
+  return nothing
+end
+
+
+function prolong2mortars!(cache, cache_gradients, u::AbstractArray{<:Any,4}, gradients,
+                          equations, mortar_l2::LobattoLegendreMortarL2, dg::DGSEM)
+  prolong2mortars!(cache, u, equations, mortar_l2, dg)
+  prolong2mortars!(cache_gradients[1], gradients[1], equations, mortar_l2, dg)
+  prolong2mortars!(cache_gradients[2], gradients[2], equations, mortar_l2, dg)
+
+  return nothing
+end
+
+
+function calc_mortar_flux!(surface_flux_values::AbstractArray{<:Any,4},
+                           nonconservative_terms::Val{false}, equations,
+                           mortar_l2::LobattoLegendreMortarL2, dg::DG, cache, cache_gradients)
+  @unpack neighbor_ids, u_lower, u_upper, orientations = cache.mortars
+  @unpack fstar_upper_threaded, fstar_lower_threaded = cache
+  gradients_x_upper = cache_gradients[1].mortars.u_upper
+  gradients_x_lower = cache_gradients[1].mortars.u_lower
+  gradients_y_upper = cache_gradients[2].mortars.u_upper
+  gradients_y_lower = cache_gradients[2].mortars.u_lower
+  gradients_upper = (gradients_x_upper, gradients_y_upper)
+  gradients_lower = (gradients_x_lower, gradients_y_lower)
+
+  Threads.@threads for mortar in eachmortar(dg, cache)
+    # Choose thread-specific pre-allocated container
+    fstar_upper = fstar_upper_threaded[Threads.threadid()]
+    fstar_lower = fstar_lower_threaded[Threads.threadid()]
+
+    # Calculate fluxes
+    orientation = orientations[mortar]
+    calc_fstar!(fstar_upper, equations, dg, u_upper, gradients_upper, mortar, orientation)
+    calc_fstar!(fstar_lower, equations, dg, u_lower, gradients_lower, mortar, orientation)
+
+    mortar_fluxes_to_elements!(surface_flux_values, equations, mortar_l2, dg, cache,
+                               mortar, fstar_upper, fstar_lower)
+  end
+
+  return nothing
+end
+
+
+@inline function calc_fstar!(destination::AbstractArray{<:Any,2}, equations, dg::DGSEM,
+                             u_interfaces, gradients_interfaces, interface, orientation)
+  @unpack surface_flux = dg
+
+  for i in eachnode(dg)
+    # Call pointwise two-point numerical flux function
+    u_ll, u_rr = get_surface_node_vars(u_interfaces, equations, dg, i, interface)
+    gradients_x_ll, gradients_x_rr = get_surface_node_vars(gradients_interfaces[1], equations, dg, i, interface)
+    gradients_y_ll, gradients_y_rr = get_surface_node_vars(gradients_interfaces[2], equations, dg, i, interface)
+    gradients_ll = (gradients_x_ll, gradients_y_ll)
+    gradients_rr = (gradients_x_rr, gradients_y_rr)
+    flux = surface_flux(u_ll, u_rr, gradients_ll, gradients_rr, orientation, equations)
+
+    # Copy flux to left and right element storage
+    set_node_vars!(destination, flux, equations, dg, i)
+  end
+
+  return nothing
+end
diff --git a/src/visualization/plot_recipes.jl b/src/visualization/plot_recipes.jl
index 566b9116f59..89cfa4a5647 100644
--- a/src/visualization/plot_recipes.jl
+++ b/src/visualization/plot_recipes.jl
@@ -249,9 +249,9 @@ end
   legend -->  :none
 
   # Set series properties
-  seriestype := :path
-  linecolor := :black
-  linewidth := 1
+  seriestype --> :path
+  linecolor --> :black
+  linewidth --> 1
 
   # Return data for plotting
   mesh_vertices_x, mesh_vertices_y
diff --git a/test/test_examples_2d_advdiff.jl b/test/test_examples_2d_advdiff.jl
new file mode 100644
index 00000000000..e29c6806718
--- /dev/null
+++ b/test/test_examples_2d_advdiff.jl
@@ -0,0 +1,25 @@
+module TestExamples2DAdvDiff
+
+using Test
+using Trixi
+
+include("test_trixi.jl")
+
+# pathof(Trixi) returns /path/to/Trixi/src/Trixi.jl, dirname gives the parent directory
+EXAMPLES_DIR = joinpath(pathof(Trixi) |> dirname |> dirname, "examples", "2d")
+
+@testset "Scalar advection-diffusion" begin
+  @testset "elixir_advdiff_basic.jl" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advdiff_basic.jl"),
+      l2   = [6.707802391539727e-5],
+      linf = [0.0005912935421852339])
+  end
+
+  @testset "elixir_advdiff_amr.jl" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advdiff_amr.jl"),
+      l2   = [0.16503674544814984],
+      linf = [0.9236893624634693])
+  end
+end
+
+end # module
diff --git a/test/test_examples_2d_heat.jl b/test/test_examples_2d_heat.jl
new file mode 100644
index 00000000000..0d5d0388377
--- /dev/null
+++ b/test/test_examples_2d_heat.jl
@@ -0,0 +1,37 @@
+module TestExamples2DHeat
+
+using Test
+using Trixi
+
+include("test_trixi.jl")
+
+# pathof(Trixi) returns /path/to/Trixi/src/Trixi.jl, dirname gives the parent directory
+EXAMPLES_DIR = joinpath(pathof(Trixi) |> dirname |> dirname, "examples", "2d")
+
+@testset "Heat equation" begin
+  @testset "elixir_heat_basic.jl" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_heat_basic.jl"),
+      l2   = [1.2700756869568872e-8],
+      linf = [6.629198634477973e-8])
+  end
+
+  @testset "elixir_heat_amr.jl" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_heat_amr.jl"),
+      l2   = [0.08891001605654104],
+      linf = [0.6575322065283009])
+  end
+
+  @testset "elixir_heat_nonperiodic.jl" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_heat_nonperiodic.jl"),
+      l2   = [2.050593060856897e-8],
+      linf = [6.631714299464696e-8])
+  end
+
+  @testset "elixir_heat_source_terms.jl" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_heat_source_terms.jl"),
+      l2   = [0.025350269628084208],
+      linf = [0.09030986815695541])
+  end
+end
+
+end # module
diff --git a/test/test_examples_2d_part1.jl b/test/test_examples_2d_part1.jl
index 029f516cfd3..d68a6beeb5e 100644
--- a/test/test_examples_2d_part1.jl
+++ b/test/test_examples_2d_part1.jl
@@ -19,6 +19,12 @@ isdir(outdir) && rm(outdir, recursive=true)
   # Linear advection
   include("test_examples_2d_advection.jl")
 
+  # Linear advection-diffusion
+  include("test_examples_2d_advdiff.jl")
+
+  # Heat equation
+  include("test_examples_2d_heat.jl")
+
   # Hyperbolic diffusion
   include("test_examples_2d_hypdiff.jl")