diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index f257da07a44..6aba503691e 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -69,6 +69,7 @@ jobs:
           - p4est_part1
           - p4est_part2
           - unstructured_dgmulti
+          - parabolic
           - paper_self_gravitating_gas_dynamics
           - misc_part1
           - misc_part2
diff --git a/examples/dgmulti_2d/elixir_advection_diffusion.jl b/examples/dgmulti_2d/elixir_advection_diffusion.jl
new file mode 100644
index 00000000000..21effda986e
--- /dev/null
+++ b/examples/dgmulti_2d/elixir_advection_diffusion.jl
@@ -0,0 +1,56 @@
+using Trixi, OrdinaryDiffEq
+
+dg = DGMulti(polydeg = 3, element_type = Tri(), approximation_type = Polynomial(),
+             surface_integral = SurfaceIntegralWeakForm(flux_lax_friedrichs),
+             volume_integral = VolumeIntegralWeakForm())
+
+equations = LinearScalarAdvectionEquation2D(1.5, 1.0)
+equations_parabolic = LaplaceDiffusion2D(5e-2, equations)
+
+initial_condition_zero(x, t, equations::LinearScalarAdvectionEquation2D) = SVector(0.0)
+initial_condition = initial_condition_zero
+
+# tag different boundary segments
+left(x, tol=50*eps()) = abs(x[1] + 1) < tol
+right(x, tol=50*eps()) = abs(x[1] - 1) < tol
+bottom(x, tol=50*eps()) = abs(x[2] + 1) < tol
+top(x, tol=50*eps()) = abs(x[2] - 1) < tol
+is_on_boundary = Dict(:left => left, :right => right, :top => top, :bottom => bottom)
+mesh = DGMultiMesh(dg, cells_per_dimension=(16, 16); is_on_boundary)
+
+# BC types
+boundary_condition_left = BoundaryConditionDirichlet((x, t, equations) -> SVector(1 + 0.1 * x[2]))
+boundary_condition_zero = BoundaryConditionDirichlet((x, t, equations) -> SVector(0.0))
+boundary_condition_neumann_zero = BoundaryConditionNeumann((x, t, equations) -> SVector(0.0))
+
+# define inviscid boundary conditions
+boundary_conditions = (; :left => boundary_condition_left,
+                         :bottom => boundary_condition_zero,
+                         :top => boundary_condition_do_nothing,
+                         :right => boundary_condition_do_nothing)
+
+# define viscous boundary conditions
+boundary_conditions_parabolic = (; :left => boundary_condition_left,
+                                   :bottom => boundary_condition_zero,
+                                   :top => boundary_condition_zero,
+                                   :right => boundary_condition_neumann_zero)
+
+semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), initial_condition, dg;
+                                             boundary_conditions=(boundary_conditions, boundary_conditions_parabolic))
+
+tspan = (0.0, 1.5)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+alive_callback = AliveCallback(alive_interval=10)
+analysis_interval = 100
+analysis_callback = AnalysisCallback(semi, interval=analysis_interval, uEltype=real(dg))
+callbacks = CallbackSet(summary_callback, alive_callback)
+
+###############################################################################
+# run the simulation
+
+time_int_tol = 1e-6
+sol = solve(ode, RDPK3SpFSAL49(), abstol=time_int_tol, reltol=time_int_tol,
+            save_everystep=false, callback=callbacks)
+summary_callback() # print the timer summary
diff --git a/examples/dgmulti_2d/elixir_advection_diffusion_periodic.jl b/examples/dgmulti_2d/elixir_advection_diffusion_periodic.jl
new file mode 100644
index 00000000000..ceeb01f71f1
--- /dev/null
+++ b/examples/dgmulti_2d/elixir_advection_diffusion_periodic.jl
@@ -0,0 +1,36 @@
+
+using Trixi, OrdinaryDiffEq
+
+dg = DGMulti(polydeg = 3, element_type = Tri(), approximation_type = Polynomial(),
+             surface_integral = SurfaceIntegralWeakForm(flux_lax_friedrichs),
+             volume_integral = VolumeIntegralWeakForm())
+
+equations = LinearScalarAdvectionEquation2D(1.0, 1.0)
+equations_parabolic = LaplaceDiffusion2D(1e-2, equations)
+
+function initial_condition_sharp_gaussian(x, t, equations::LinearScalarAdvectionEquation2D)
+  return SVector(exp(-100 * (x[1]^2 + x[2]^2)))
+end
+initial_condition = initial_condition_sharp_gaussian
+
+mesh = DGMultiMesh(dg, cells_per_dimension = (16, 16), periodicity=true)
+semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic),
+                                             initial_condition, dg)
+
+tspan = (0.0, 2.0)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+alive_callback = AliveCallback(alive_interval=10)
+analysis_interval = 100
+analysis_callback = AnalysisCallback(semi, interval=analysis_interval, uEltype=real(dg))
+callbacks = CallbackSet(summary_callback, alive_callback, analysis_callback)
+
+###############################################################################
+# run the simulation
+
+time_int_tol = 1e-6
+sol = solve(ode, RDPK3SpFSAL49(), abstol=time_int_tol, reltol=time_int_tol,
+            save_everystep=false, callback=callbacks)
+
+summary_callback() # print the timer summary
diff --git a/src/Trixi.jl b/src/Trixi.jl
index 31c8a73662c..5ee6e7afb95 100644
--- a/src/Trixi.jl
+++ b/src/Trixi.jl
@@ -26,7 +26,8 @@ using SparseArrays: AbstractSparseMatrix, AbstractSparseMatrixCSC, sparse, dropt
 using Reexport: @reexport
 
 using SciMLBase: CallbackSet, DiscreteCallback,
-                 ODEProblem, ODESolution, ODEFunction
+                 ODEProblem, ODESolution, ODEFunction,
+                 SplitODEProblem
 import SciMLBase: get_du, get_tmp_cache, u_modified!,
                   AbstractODEIntegrator, init, step!, check_error,
                   get_proposed_dt, set_proposed_dt!,
@@ -106,6 +107,7 @@ include("meshes/meshes.jl")
 include("solvers/solvers.jl")
 include("semidiscretization/semidiscretization.jl")
 include("semidiscretization/semidiscretization_hyperbolic.jl")
+include("semidiscretization/semidiscretization_hyperbolic_parabolic.jl")
 include("semidiscretization/semidiscretization_euler_acoustics.jl")
 include("callbacks_step/callbacks_step.jl")
 include("callbacks_stage/callbacks_stage.jl")
@@ -128,6 +130,7 @@ export AcousticPerturbationEquations2D,
        HyperbolicDiffusionEquations1D, HyperbolicDiffusionEquations2D, HyperbolicDiffusionEquations3D,
        LinearScalarAdvectionEquation1D, LinearScalarAdvectionEquation2D, LinearScalarAdvectionEquation3D,
        InviscidBurgersEquation1D,
+       LaplaceDiffusion2D,
        LatticeBoltzmannEquations2D, LatticeBoltzmannEquations3D,
        ShallowWaterEquations1D, ShallowWaterEquations2D
 
@@ -151,8 +154,10 @@ export initial_condition_constant,
        initial_condition_density_wave,
        initial_condition_weak_blast_wave
 
-export boundary_condition_periodic,
+export boundary_condition_do_nothing,
+       boundary_condition_periodic,
        BoundaryConditionDirichlet,
+       BoundaryConditionNeumann,
        boundary_condition_noslip_wall,
        boundary_condition_slip_wall,
        boundary_condition_wall
@@ -185,6 +190,8 @@ export nelements, nnodes, nvariables,
 
 export SemidiscretizationHyperbolic, semidiscretize, compute_coefficients, integrate
 
+export SemidiscretizationHyperbolicParabolic
+
 export SemidiscretizationEulerAcoustics
 
 export SemidiscretizationEulerGravity, ParametersEulerGravity,
diff --git a/src/equations/equations.jl b/src/equations/equations.jl
index dc3bfffba30..7cc6768d116 100644
--- a/src/equations/equations.jl
+++ b/src/equations/equations.jl
@@ -175,6 +175,45 @@ end
   return flux
 end
 
+# operator types used for dispatch on parabolic boundary fluxes
+struct Gradient end
+struct Divergence end
+
+# Dirichlet-type boundary condition for use with parabolic equations.
+@inline function (boundary_condition::BoundaryConditionDirichlet)(u_inner, normal::AbstractVector,
+                                                                  x, t, operator_type::Divergence,
+                                                                  equations)
+  return u_inner
+end
+
+
+"""
+    BoundaryConditionNeumann(boundary_normal_flux_function)
+
+Similar to `BoundaryConditionDirichlet`, but creates a Neumann boundary condition for parabolic
+equations that uses the function `boundary_normal_flux_function` to specify the values of the normal
+flux at the boundary.
+The passed boundary value function will be called with the same arguments as an initial condition function is called, i.e., as
+```julia
+boundary_normal_flux_function(x, t, equations)
+```
+where `x` specifies the coordinates, `t` is the current time, and `equation` is the corresponding system of equations.
+"""
+struct BoundaryConditionNeumann{B}
+  boundary_normal_flux_function::B
+end
+
+# specify default behavior for a Neumann BC is to evaluate the flux at the inner state
+@inline function (boundary_condition::BoundaryConditionNeumann)(flux_inner, normal::AbstractVector,
+                                                                x, t, operator_type::Gradient,
+                                                                equations)
+  return flux_inner
+end
+
+# Imposing no boundary condition just evaluates the flux at the inner state.
+@inline function boundary_condition_do_nothing(inner_flux_or_state, other_args...)
+  return inner_flux_or_state
+end
 
 # set sensible default values that may be overwritten by specific equations
 """
@@ -329,5 +368,10 @@ include("lattice_boltzmann_3d.jl")
 abstract type AbstractAcousticPerturbationEquations{NDIMS, NVARS} <: AbstractEquations{NDIMS, NVARS} end
 include("acoustic_perturbation_2d.jl")
 
+abstract type AbstractEquationsParabolic{NDIMS, NVARS} <: AbstractEquations{NDIMS, NVARS} end
+
+# Linear scalar diffusion for use in linear scalar advection-diffusion problems
+abstract type AbstractLaplaceDiffusionEquations{NDIMS, NVARS} <: AbstractEquationsParabolic{NDIMS, NVARS} end
+include("laplace_diffusion_2d.jl")
 
 end # @muladd
diff --git a/src/equations/laplace_diffusion_2d.jl b/src/equations/laplace_diffusion_2d.jl
new file mode 100644
index 00000000000..774d93cf7b3
--- /dev/null
+++ b/src/equations/laplace_diffusion_2d.jl
@@ -0,0 +1,32 @@
+@doc raw"""
+    LaplaceDiffusion2D(diffusivity, equations)
+
+`LaplaceDiffusion2D` represents a scalar diffusion term ``\nabla \cdot (\kappa\nabla u))``
+with diffusivity ``\kappa`` applied to each solution component defined by `equations`.
+"""
+struct LaplaceDiffusion2D{E, N, T} <: AbstractLaplaceDiffusionEquations{2, N}
+  diffusivity::T
+  equations::E
+end
+
+LaplaceDiffusion2D(diffusivity, equations) =
+  LaplaceDiffusion2D{typeof(equations), nvariables(equations), typeof(diffusivity)}(diffusivity, equations)
+
+# no orientation specified since the flux is vector-valued
+function flux(u, grad_u, equations::LaplaceDiffusion2D)
+  dudx, dudy = grad_u
+  return SVector(equations.diffusivity * dudx, equations.diffusivity * dudy)
+end
+
+# Dirichlet-type boundary condition for use with a parabolic solver in weak form
+@inline function (boundary_condition::BoundaryConditionDirichlet)(u_inner, normal::AbstractVector,
+                                                                  x, t, operator_type::Gradient,
+                                                                  equations::LaplaceDiffusion2D)
+  return boundary_condition.boundary_value_function(x, t, equations)
+end
+
+@inline function (boundary_condition::BoundaryConditionNeumann)(flux_inner, normal::AbstractVector,
+                                                                x, t, operator_type::Divergence,
+                                                                equations::LaplaceDiffusion2D)
+  return boundary_condition.boundary_normal_flux_function(x, t, equations)
+end
\ No newline at end of file
diff --git a/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl b/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl
new file mode 100644
index 00000000000..ea14b05e478
--- /dev/null
+++ b/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl
@@ -0,0 +1,256 @@
+# By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
+# Since these FMAs can increase the performance of many numerical algorithms,
+# we need to opt-in explicitly.
+# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
+@muladd begin
+
+
+"""
+    SemidiscretizationHyperbolicParabolic
+
+A struct containing everything needed to describe a spatial semidiscretization
+of a mixed hyperbolic-parabolic conservation law.
+"""
+struct SemidiscretizationHyperbolicParabolic{Mesh, Equations, EquationsParabolic, InitialCondition,
+                                             BoundaryConditions, BoundaryConditionsParabolic,
+                                             SourceTerms, Solver, Cache, CacheParabolic} <: AbstractSemidiscretization
+
+  mesh::Mesh
+
+  equations::Equations
+  equations_parabolic::EquationsParabolic
+
+  # 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
+  boundary_conditions_parabolic::BoundaryConditionsParabolic
+
+  source_terms::SourceTerms
+
+  solver::Solver
+  # TODO: introduce `solver_parabolic` for future specialization.
+
+  cache::Cache
+  cache_parabolic::CacheParabolic
+
+  performance_counter::PerformanceCounter
+
+  function SemidiscretizationHyperbolicParabolic{Mesh, Equations, EquationsParabolic, InitialCondition, BoundaryConditions, BoundaryConditionsParabolic, SourceTerms, Solver, Cache, CacheParabolic}(
+      mesh::Mesh, equations::Equations, equations_parabolic::EquationsParabolic, initial_condition::InitialCondition,
+      boundary_conditions::BoundaryConditions, boundary_conditions_parabolic::BoundaryConditionsParabolic,
+      source_terms::SourceTerms, solver::Solver, cache::Cache, cache_parabolic::CacheParabolic) where {Mesh, Equations, EquationsParabolic, InitialCondition, BoundaryConditions, BoundaryConditionsParabolic, SourceTerms, Solver, Cache, CacheParabolic}
+    @assert ndims(mesh) == ndims(equations)
+
+    # Todo: assert nvariables(equations)==nvariables(equations_parabolic)
+
+    performance_counter = PerformanceCounter()
+
+    new(mesh, equations, equations_parabolic, initial_condition,
+        boundary_conditions, boundary_conditions_parabolic,
+        source_terms, solver, cache, cache_parabolic, performance_counter)
+  end
+end
+
+"""
+    SemidiscretizationHyperbolicParabolic(mesh, both_equations, initial_condition, solver;
+                                          source_terms=nothing,
+                                          both_boundary_conditions=(boundary_condition_periodic, boundary_condition_periodic),
+                                          RealT=real(solver),
+                                          uEltype=RealT,
+                                          both_initial_caches=(NamedTuple(), NamedTuple()))
+
+Construct a semidiscretization of a hyperbolic-parabolic PDE.
+"""
+function SemidiscretizationHyperbolicParabolic(mesh, equations::Tuple,
+                                               initial_condition, solver;
+                                               source_terms=nothing,
+                                               boundary_conditions=(boundary_condition_periodic, boundary_condition_periodic),
+                                               # `RealT` is used as real type for node locations etc.
+                                               # while `uEltype` is used as element type of solutions etc.
+                                               RealT=real(solver), uEltype=RealT,
+                                               initial_caches=(NamedTuple(), NamedTuple()))
+
+  equations_hyperbolic, equations_parabolic = equations
+  boundary_conditions_hyperbolic, boundary_conditions_parabolic = boundary_conditions
+  initial_hyperbolic_cache, initial_cache_parabolic = initial_caches
+
+  return SemidiscretizationHyperbolicParabolic(mesh, equations_hyperbolic, equations_parabolic,
+                                               initial_condition, solver; source_terms,
+                                               boundary_conditions=boundary_conditions_hyperbolic,
+                                               boundary_conditions_parabolic=boundary_conditions_parabolic,
+                                               RealT, uEltype, initial_cache=initial_hyperbolic_cache,
+                                               initial_cache_parabolic=initial_cache_parabolic)
+end
+
+function SemidiscretizationHyperbolicParabolic(mesh, equations, equations_parabolic,
+                                               initial_condition, solver;
+                                               source_terms=nothing,
+                                               boundary_conditions=boundary_condition_periodic,
+                                               boundary_conditions_parabolic=boundary_condition_periodic,
+                                               # `RealT` is used as real type for node locations etc.
+                                               # while `uEltype` is used as element type of solutions etc.
+                                               RealT=real(solver), uEltype=RealT,
+                                               initial_cache=NamedTuple(),
+                                               initial_cache_parabolic=NamedTuple())
+
+  cache = (; create_cache(mesh, equations, solver, RealT, uEltype)..., initial_cache...)
+  _boundary_conditions = digest_boundary_conditions(boundary_conditions, mesh, solver, cache)
+  _boundary_conditions_parabolic = digest_boundary_conditions(boundary_conditions_parabolic, mesh, solver, cache)
+
+  cache_parabolic = (; create_cache_parabolic(mesh, equations_parabolic, solver, RealT, uEltype)...,
+                                              initial_cache_parabolic...)
+
+  SemidiscretizationHyperbolicParabolic{typeof(mesh), typeof(equations), typeof(equations_parabolic),
+                                        typeof(initial_condition), typeof(_boundary_conditions), typeof(_boundary_conditions_parabolic),
+                                        typeof(source_terms), typeof(solver), typeof(cache), typeof(cache_parabolic)}(
+    mesh, equations, equations_parabolic, initial_condition,
+    _boundary_conditions, _boundary_conditions_parabolic, source_terms,
+    solver, cache, cache_parabolic)
+end
+
+
+# Create a new semidiscretization but change some parameters compared to the input.
+# `Base.similar` follows a related concept but would require us to `copy` the `mesh`,
+# which would impact the performance. Instead, `SciMLBase.remake` has exactly the
+# semantics we want to use here. In particular, it allows us to re-use mutable parts,
+# e.g. `remake(semi).mesh === semi.mesh`.
+function remake(semi::SemidiscretizationHyperbolicParabolic; uEltype=real(semi.solver),
+                                                             mesh=semi.mesh,
+                                                             equations=semi.equations,
+                                                             equations_parabolic=semi.equations_parabolic,
+                                                             initial_condition=semi.initial_condition,
+                                                             solver=semi.solver,
+                                                             source_terms=semi.source_terms,
+                                                             boundary_conditions=semi.boundary_conditions,
+                                                             boundary_conditions_parabolic=semi.boundary_conditions_parabolic
+                                                             )
+  # TODO: Which parts do we want to `remake`? At least the solver needs some
+  #       special care if shock-capturing volume integrals are used (because of
+  #       the indicators and their own caches...).
+  SemidiscretizationHyperbolicParabolic(
+    mesh, equations, equations_parabolic, initial_condition, solver; source_terms, boundary_conditions, boundary_conditions_parabolic, uEltype)
+end
+
+function Base.show(io::IO, semi::SemidiscretizationHyperbolicParabolic)
+  @nospecialize semi # reduce precompilation time
+
+  print(io, "SemidiscretizationHyperbolicParabolic(")
+  print(io,       semi.mesh)
+  print(io, ", ", semi.equations)
+  print(io, ", ", semi.equations_parabolic)
+  print(io, ", ", semi.initial_condition)
+  print(io, ", ", semi.boundary_conditions)
+  print(io, ", ", semi.boundary_conditions_parabolic)
+  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::SemidiscretizationHyperbolicParabolic)
+  @nospecialize semi # reduce precompilation time
+
+  if get(io, :compact, false)
+    show(io, semi)
+  else
+    summary_header(io, "SemidiscretizationHyperbolicParabolic")
+    summary_line(io, "#spatial dimensions", ndims(semi.equations))
+    summary_line(io, "mesh", semi.mesh)
+    summary_line(io, "hyperbolic equations", semi.equations |> typeof |> nameof)
+    summary_line(io, "parabolic equations", semi.equations_parabolic |> typeof |> nameof)
+    summary_line(io, "initial condition", semi.initial_condition)
+
+    # print_boundary_conditions(io, semi)
+
+    summary_line(io, "source terms", semi.source_terms)
+    summary_line(io, "solver", semi.solver |> typeof |> nameof)
+    summary_line(io, "total #DOFs", ndofs(semi))
+    summary_footer(io)
+  end
+end
+
+@inline Base.ndims(semi::SemidiscretizationHyperbolicParabolic) = ndims(semi.mesh)
+
+@inline nvariables(semi::SemidiscretizationHyperbolicParabolic) = nvariables(semi.equations)
+
+@inline Base.real(semi::SemidiscretizationHyperbolicParabolic) = real(semi.solver)
+
+# retain dispatch on hyperbolic equations only
+@inline function mesh_equations_solver_cache(semi::SemidiscretizationHyperbolicParabolic)
+  @unpack mesh, equations, solver, cache = semi
+  return mesh, equations, solver, cache
+end
+
+
+function calc_error_norms(func, u_ode, t, analyzer, semi::SemidiscretizationHyperbolicParabolic, 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 compute_coefficients(t, semi::SemidiscretizationHyperbolicParabolic)
+  # Call `compute_coefficients` in `src/semidiscretization/semidiscretization.jl`
+  compute_coefficients(semi.initial_condition, t, semi)
+end
+
+function compute_coefficients!(u_ode, t, semi::SemidiscretizationHyperbolicParabolic)
+  compute_coefficients!(u_ode, semi.initial_condition, t, semi)
+end
+
+"""
+    semidiscretize(semi::SemidiscretizationHyperbolicParabolic, tspan)
+
+Wrap the semidiscretization `semi` as a Split ODE problem in the time interval `tspan`
+that can be passed to `solve` from the [SciML ecosystem](https://diffeq.sciml.ai/latest/).
+"""
+function semidiscretize(semi::SemidiscretizationHyperbolicParabolic, tspan)
+  u0_ode = compute_coefficients(first(tspan), semi)
+  # TODO: MPI, do we want to synchonize loading and print debug statements, e.g. using
+  #       mpi_isparallel() && MPI.Barrier(mpi_comm())
+  #       See https://github.com/trixi-framework/Trixi.jl/issues/328
+  return SplitODEProblem(rhs!, rhs_parabolic!, u0_ode, tspan, semi)
+end
+
+function rhs!(du_ode, u_ode, semi::SemidiscretizationHyperbolicParabolic, 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)
+
+  # TODO: Taal decide, do we need to pass the mesh?
+  time_start = time_ns()
+  @trixi_timeit timer() "rhs!" rhs!(du, u, t, mesh, equations, initial_condition,
+                                    boundary_conditions, source_terms, solver, cache)
+  runtime = time_ns() - time_start
+  put!(semi.performance_counter, runtime)
+
+  return nothing
+end
+
+function rhs_parabolic!(du_ode, u_ode, semi::SemidiscretizationHyperbolicParabolic, t)
+  @unpack mesh, equations_parabolic, initial_condition, boundary_conditions_parabolic, source_terms, solver, cache, cache_parabolic = semi
+
+  u  = wrap_array(u_ode,  mesh, equations_parabolic, solver, cache_parabolic)
+  du = wrap_array(du_ode, mesh, equations_parabolic, solver, cache_parabolic)
+
+  # TODO: Taal decide, do we need to pass the mesh?
+  time_start = time_ns()
+  @trixi_timeit timer() "parabolic rhs!" rhs_parabolic!(du, u, t, mesh, equations_parabolic, initial_condition,
+                                                        boundary_conditions_parabolic, source_terms,
+                                                        solver, cache, cache_parabolic)
+  runtime = time_ns() - time_start
+  put!(semi.performance_counter, runtime)
+
+  return nothing
+end
+
+
+end # @muladd
diff --git a/src/solvers/dgmulti.jl b/src/solvers/dgmulti.jl
index d5c8f12d395..318a11b678e 100644
--- a/src/solvers/dgmulti.jl
+++ b/src/solvers/dgmulti.jl
@@ -9,3 +9,6 @@ include("dgmulti/sbp.jl")
 
 # specialization of DGMulti to specific equations
 include("dgmulti/flux_differencing_compressible_euler.jl")
+
+# parabolic terms for DGMulti solvers
+include("dgmulti/dg_parabolic.jl")
\ No newline at end of file
diff --git a/src/solvers/dgmulti/dg.jl b/src/solvers/dgmulti/dg.jl
index e4d25920eac..6f2a835dbd4 100644
--- a/src/solvers/dgmulti/dg.jl
+++ b/src/solvers/dgmulti/dg.jl
@@ -40,24 +40,24 @@ end
 @inline eachdim(mesh) = Base.OneTo(ndims(mesh))
 
 # iteration over all elements in a mesh
-@inline ndofs(mesh::DGMultiMesh, dg::DGMulti, cache) = dg.basis.Np * mesh.md.num_elements
-@inline eachelement(mesh::DGMultiMesh, dg::DGMulti, cache) = Base.OneTo(mesh.md.num_elements)
+@inline ndofs(mesh::DGMultiMesh, dg::DGMulti, other_args...) = dg.basis.Np * mesh.md.num_elements
+@inline eachelement(mesh::DGMultiMesh, dg::DGMulti, other_args...) = Base.OneTo(mesh.md.num_elements)
 
 # iteration over quantities in a single element
 @inline nnodes(basis::RefElemData) = basis.Np
-@inline each_face_node(mesh::DGMultiMesh, dg::DGMulti, cache) = Base.OneTo(dg.basis.Nfq)
-@inline each_quad_node(mesh::DGMultiMesh, dg::DGMulti, cache) = Base.OneTo(dg.basis.Nq)
+@inline each_face_node(mesh::DGMultiMesh, dg::DGMulti, other_args...) = Base.OneTo(dg.basis.Nfq)
+@inline each_quad_node(mesh::DGMultiMesh, dg::DGMulti, other_args...) = Base.OneTo(dg.basis.Nq)
 
 # iteration over quantities over the entire mesh (dofs, quad nodes, face nodes).
-@inline each_dof_global(mesh::DGMultiMesh, dg::DGMulti, cache) = Base.OneTo(ndofs(mesh, dg, cache))
-@inline each_quad_node_global(mesh::DGMultiMesh, dg::DGMulti, cache) = Base.OneTo(dg.basis.Nq * mesh.md.num_elements)
-@inline each_face_node_global(mesh::DGMultiMesh, dg::DGMulti, cache) = Base.OneTo(dg.basis.Nfq * mesh.md.num_elements)
+@inline each_dof_global(mesh::DGMultiMesh, dg::DGMulti, other_args...) = Base.OneTo(ndofs(mesh, dg, other_args...))
+@inline each_quad_node_global(mesh::DGMultiMesh, dg::DGMulti, other_args...) = Base.OneTo(dg.basis.Nq * mesh.md.num_elements)
+@inline each_face_node_global(mesh::DGMultiMesh, dg::DGMulti, other_args...) = Base.OneTo(dg.basis.Nfq * mesh.md.num_elements)
 
 # interface with semidiscretization_hyperbolic
 wrap_array(u_ode, mesh::DGMultiMesh, equations, dg::DGMulti, cache) = u_ode
 wrap_array_native(u_ode, mesh::DGMultiMesh, equations, dg::DGMulti, cache) = u_ode
-function digest_boundary_conditions(boundary_conditions::NamedTuple{Keys,ValueTypes}, mesh::DGMultiMesh,
-                                    dg::DGMulti, cache) where {Keys,ValueTypes<:NTuple{N,Any}} where {N}
+function digest_boundary_conditions(boundary_conditions::NamedTuple{Keys, ValueTypes}, mesh::DGMultiMesh,
+                                    dg::DGMulti, cache) where {Keys, ValueTypes<:NTuple{N, Any}} where {N}
   return boundary_conditions
 end
 
@@ -67,7 +67,7 @@ function allocate_nested_array(uEltype, nvars, array_dimensions, dg)
   return StructArray{SVector{nvars, uEltype}}(ntuple(_->zeros(uEltype, array_dimensions...), nvars))
 end
 
-function reset_du!(du, dg::DGMulti, cache)
+function reset_du!(du, dg::DGMulti, other_args...)
   @threaded for i in eachindex(du)
       du[i] = zero(eltype(du))
   end
@@ -208,13 +208,9 @@ function calc_interface_flux!(cache, surface_integral::SurfaceIntegralWeakForm,
     # inner (idM -> minus) and outer (idP -> plus) indices
     idM, idP = mapM[face_node_index], mapP[face_node_index]
     uM = u_face_values[idM]
-
-    # compute flux if node is not a boundary node
-    if idM != idP
-      uP = u_face_values[idP]
-      normal = SVector{NDIMS}(getindex.(nxyzJ, idM)) / Jf[idM]
-      flux_face_values[idM] = surface_flux(uM, uP, normal, equations) * Jf[idM]
-    end
+    uP = u_face_values[idP]
+    normal = SVector{NDIMS}(getindex.(nxyzJ, idM)) / Jf[idM]
+    flux_face_values[idM] = surface_flux(uM, uP, normal, equations) * Jf[idM]
   end
 end
 
@@ -353,9 +349,10 @@ end
 
 
 # Todo: DGMulti. Specialize for modal DG on curved meshes using WADG
-function invert_jacobian!(du, mesh::DGMultiMesh, equations, dg::DGMulti, cache)
+# inverts Jacobian and scales by -1.0
+function invert_jacobian!(du, mesh::DGMultiMesh, equations, dg::DGMulti, cache; scaling=-1)
   @threaded for i in each_dof_global(mesh, dg, cache)
-    du[i] *= -cache.invJ[i]
+    du[i] *= scaling * cache.invJ[i]
   end
 end
 
diff --git a/src/solvers/dgmulti/dg_parabolic.jl b/src/solvers/dgmulti/dg_parabolic.jl
new file mode 100644
index 00000000000..f4c4799e021
--- /dev/null
+++ b/src/solvers/dgmulti/dg_parabolic.jl
@@ -0,0 +1,275 @@
+function create_cache_parabolic(mesh::DGMultiMesh, equations::AbstractEquationsParabolic,
+                                dg::DGMulti, RealT, uEltype)
+  nvars = nvariables(equations)
+
+  @unpack M, Drst = dg.basis
+  weak_differentiation_matrices = map(A -> -M \ (A' * M), Drst)
+
+  # u_transformed stores "transformed" variables for computing the gradient
+  @unpack md = mesh
+  u_transformed = allocate_nested_array(uEltype, nvars, size(md.x), dg)
+  u_grad = ntuple(_ -> similar(u_transformed), ndims(mesh))
+  viscous_flux = similar.(u_grad)
+
+  u_face_values = allocate_nested_array(uEltype, nvars, size(md.xf), dg)
+  scalar_flux_face_values = similar(u_face_values)
+  grad_u_face_values = ntuple(_ -> similar(u_face_values), ndims(mesh))
+
+  local_u_values_threaded = [similar(u_transformed, dg.basis.Nq) for _ in 1:Threads.nthreads()]
+  local_viscous_flux_threaded = [ntuple(_ -> similar(u_transformed, dg.basis.Nq), ndims(mesh)) for _ in 1:Threads.nthreads()]
+  local_flux_face_values_threaded = [similar(scalar_flux_face_values[:, 1]) for _ in 1:Threads.nthreads()]
+
+  return (; u_transformed, u_grad, viscous_flux,
+            weak_differentiation_matrices,
+            u_face_values, grad_u_face_values, scalar_flux_face_values,
+            local_u_values_threaded, local_viscous_flux_threaded, local_flux_face_values_threaded)
+end
+
+# Transform solution variables prior to taking the gradient
+# (e.g., conservative to primitive variables). Defaults to doing nothing.
+# TODO: can we avoid copying data?
+function transform_variables!(u_transformed, u, equations)
+  @threaded for i in eachindex(u)
+    u_transformed[i] = u[i]
+  end
+end
+
+# interpolates from solution coefficients to face quadrature points
+function prolong2interfaces!(u_face_values, u, mesh::DGMultiMesh, equations::AbstractEquationsParabolic,
+                             surface_integral, dg::DGMulti, cache)
+  apply_to_each_field(mul_by!(dg.basis.Vf), u_face_values, u)
+end
+
+function calc_gradient_surface_integral(u_grad, u, scalar_flux_face_values,
+                                        mesh, equations::AbstractEquationsParabolic,
+                                        dg::DGMulti, cache, cache_parabolic)
+  @unpack local_flux_face_values_threaded = cache_parabolic
+  @threaded for e in eachelement(mesh, dg)
+    local_flux_values = local_flux_face_values_threaded[Threads.threadid()]
+    for dim in eachdim(mesh)
+      for i in eachindex(local_flux_values)
+        # compute flux * (nx, ny, nz)
+        local_flux_values[i] = scalar_flux_face_values[i, e] * mesh.md.nxyzJ[dim][i, e]
+      end
+      apply_to_each_field(mul_by_accum!(dg.basis.LIFT), view(u_grad[dim], :, e), local_flux_values)
+    end
+  end
+end
+
+function calc_gradient!(u_grad, u::StructArray, t, mesh::DGMultiMesh,
+                        equations::AbstractEquationsParabolic,
+                        boundary_conditions, dg::DGMulti, cache, cache_parabolic)
+
+  @unpack weak_differentiation_matrices = cache_parabolic
+
+  for dim in eachindex(u_grad)
+    reset_du!(u_grad[dim], dg)
+  end
+
+  # compute volume contributions to gradients
+  @threaded for e in eachelement(mesh, dg)
+    for i in eachdim(mesh), j in eachdim(mesh)
+      dxidxhatj = mesh.md.rstxyzJ[i, j][1, e] # TODO: assumes mesh is affine
+      apply_to_each_field(mul_by_accum!(weak_differentiation_matrices[j], dxidxhatj),
+                          view(u_grad[i], :, e), view(u, :, e))
+    end
+  end
+
+  @unpack u_face_values = cache_parabolic
+  prolong2interfaces!(u_face_values, u, mesh, equations, dg.surface_integral, dg, cache)
+
+  # compute fluxes at interfaces
+  @unpack scalar_flux_face_values = cache_parabolic
+  @unpack mapM, mapP, Jf = mesh.md
+  @threaded for face_node_index in each_face_node_global(mesh, dg)
+    idM, idP = mapM[face_node_index], mapP[face_node_index]
+    uM = u_face_values[idM]
+    uP = u_face_values[idP]
+    scalar_flux_face_values[idM] = 0.5 * (uP + uM) # TODO: use strong/weak formulation for curved meshes?
+  end
+
+  calc_boundary_flux!(scalar_flux_face_values, nothing, t, Gradient(), boundary_conditions,
+                      mesh, equations, dg, cache, cache_parabolic)
+
+  # compute surface contributions
+  calc_gradient_surface_integral(u_grad, u, scalar_flux_face_values,
+                                 mesh, equations, dg, cache, cache_parabolic)
+
+  for dim in eachdim(mesh)
+    invert_jacobian!(u_grad[dim], mesh, equations, dg, cache; scaling=1.0)
+  end
+end
+
+# do nothing for periodic domains
+function calc_boundary_flux!(flux, u, t, operator_type, ::BoundaryConditionPeriodic,
+                             mesh, equations::AbstractEquationsParabolic, dg::DGMulti,
+                             cache, cache_parabolic)
+  return nothing
+end
+
+# "lispy tuple programming" instead of for loop for type stability
+function calc_boundary_flux!(flux, u, t, operator_type, boundary_conditions,
+                             mesh, equations, dg::DGMulti, cache, cache_parabolic)
+
+  # peel off first boundary condition
+  calc_single_boundary_flux!(flux, u, t, operator_type, first(boundary_conditions), first(keys(boundary_conditions)),
+                             mesh, equations, dg, cache, cache_parabolic)
+
+  # recurse on the remainder of the boundary conditions
+  calc_boundary_flux!(flux, u, t, operator_type, Base.tail(boundary_conditions),
+                      mesh, equations, dg, cache, cache_parabolic)
+end
+
+# terminate recursion
+calc_boundary_flux!(flux, u, t, operator_type, boundary_conditions::NamedTuple{(),Tuple{}},
+                    mesh, equations, dg::DGMulti, cache, cache_parabolic) = nothing
+
+# TODO: DGMulti. Decide if we want to use the input `u_face_values` (currently unused)
+function calc_single_boundary_flux!(flux_face_values, u_face_values, t,
+                                    operator_type, boundary_condition, boundary_key,
+                                    mesh, equations, dg::DGMulti{NDIMS}, cache, cache_parabolic) where {NDIMS}
+  rd = dg.basis
+  md = mesh.md
+
+  num_pts_per_face = rd.Nfq ÷ rd.Nfaces
+  @unpack xyzf, nxyzJ, Jf = md
+  for f in mesh.boundary_faces[boundary_key]
+    for i in Base.OneTo(num_pts_per_face)
+
+      # reverse engineer element + face node indices (avoids reshaping arrays)
+      e = ((f-1) ÷ rd.Nfaces) + 1
+      fid = i + ((f-1) % rd.Nfaces) * num_pts_per_face
+
+      face_normal = SVector{NDIMS}(getindex.(nxyzJ, fid, e)) / Jf[fid,e]
+      face_coordinates = SVector{NDIMS}(getindex.(xyzf, fid, e))
+
+      # for both the gradient and the divergence, the boundary flux is scalar valued.
+      # for the gradient, it is the solution; for divergence, it is the normal flux.
+      u_boundary = boundary_condition(flux_face_values[fid,e],
+                                      face_normal, face_coordinates, t,
+                                      operator_type, equations)
+      flux_face_values[fid,e] = u_boundary
+    end
+  end
+  return nothing
+end
+
+function calc_viscous_fluxes!(viscous_flux, u, u_grad, mesh::DGMultiMesh,
+                              equations::AbstractEquationsParabolic,
+                              dg::DGMulti, cache, cache_parabolic)
+
+  for dim in eachdim(mesh)
+    reset_du!(viscous_flux[dim], dg)
+  end
+
+  @unpack local_viscous_flux_threaded, local_u_values_threaded = cache_parabolic
+
+  @threaded for e in eachelement(mesh, dg)
+
+    # reset local storage for each element
+    local_viscous_flux = local_viscous_flux_threaded[Threads.threadid()]
+    local_u_values = local_u_values_threaded[Threads.threadid()]
+    fill!(local_u_values, zero(eltype(local_u_values)))
+    for dim in eachdim(mesh)
+      fill!(local_viscous_flux[dim], zero(eltype(local_viscous_flux[dim])))
+    end
+
+    # interpolate u and gradient to quadrature points, store in `local_viscous_flux`
+    apply_to_each_field(mul_by!(dg.basis.Vq), local_u_values, view(u, :, e)) # TODO: can we avoid this when we don't need it?
+    for dim in eachdim(mesh)
+      apply_to_each_field(mul_by!(dg.basis.Vq), local_viscous_flux[dim], view(u_grad[dim], :, e))
+    end
+
+    # compute viscous flux at quad points
+    for i in eachindex(local_u_values)
+      u_i = local_u_values[i]
+      u_grad_i = getindex.(local_viscous_flux, i) # TODO: check if this allocates. Shouldn't for tuples or SVector...
+      viscous_flux_i = flux(u_i, u_grad_i, equations)
+      setindex!.(local_viscous_flux, viscous_flux_i, i)
+    end
+
+    # project back to the DG approximation space
+    for dim in eachdim(mesh)
+      apply_to_each_field(mul_by!(dg.basis.Pq), view(viscous_flux[dim], :, e), local_viscous_flux[dim])
+    end
+  end
+end
+
+function calc_divergence!(du, u::StructArray, t, viscous_flux, mesh::DGMultiMesh,
+                          equations::AbstractEquationsParabolic,
+                          boundary_conditions, dg::DGMulti, cache, cache_parabolic)
+
+  @unpack weak_differentiation_matrices = cache_parabolic
+
+  reset_du!(du, dg)
+
+  # compute volume contributions to divergence
+  @threaded for e in eachelement(mesh, dg)
+    for i in eachdim(mesh), j in eachdim(mesh)
+      dxidxhatj = mesh.md.rstxyzJ[i, j][1, e] # assumes mesh is affine
+      apply_to_each_field(mul_by_accum!(weak_differentiation_matrices[j], dxidxhatj),
+                                view(du, :, e), view(viscous_flux[i], :, e))
+    end
+  end
+
+  # interpolates from solution coefficients to face quadrature points
+  viscous_flux_face_values = cache_parabolic.grad_u_face_values
+  for dim in eachdim(mesh)
+    prolong2interfaces!(viscous_flux_face_values[dim], viscous_flux[dim], mesh, equations,
+                        dg.surface_integral, dg, cache)
+  end
+
+  # compute fluxes at interfaces
+  @unpack scalar_flux_face_values = cache_parabolic
+  @unpack mapM, mapP, nxyzJ = mesh.md
+  @threaded for face_node_index in each_face_node_global(mesh, dg, cache, cache_parabolic)
+    idM, idP = mapM[face_node_index], mapP[face_node_index]
+
+    # compute f(u, ∇u) ⋅ n
+    flux_face_value = zero(eltype(scalar_flux_face_values))
+    for dim in eachdim(mesh)
+      uM = viscous_flux_face_values[dim][idM]
+      uP = viscous_flux_face_values[dim][idP]
+      # TODO: use strong/weak formulation?
+      flux_face_value = flux_face_value + 0.5 * (uP + uM) * nxyzJ[dim][face_node_index]
+    end
+    scalar_flux_face_values[idM] = flux_face_value
+  end
+
+  # TODO: decide what to pass in
+  calc_boundary_flux!(scalar_flux_face_values, nothing, t, Divergence(),
+                      boundary_conditions, mesh, equations, dg, cache, cache_parabolic)
+
+  # surface contributions
+  apply_to_each_field(mul_by_accum!(dg.basis.LIFT), du, scalar_flux_face_values)
+
+  invert_jacobian!(du, mesh, equations, dg, cache; scaling=1.0)
+end
+
+# assumptions: parabolic terms are of the form div(f(u, grad(u))) and
+# will be discretized first order form as follows:
+#               1. compute grad(u)
+#               2. compute f(u, grad(u))
+#               3. compute div(u)
+# boundary conditions will be applied to both grad(u) and div(u).
+function rhs_parabolic!(du, u, t, mesh::DGMultiMesh, equations_parabolic::AbstractEquationsParabolic,
+                        initial_condition, boundary_conditions, source_terms,
+                        dg::DGMulti, cache, cache_parabolic)
+
+  reset_du!(du, dg)
+
+  @unpack u_transformed, u_grad, viscous_flux = cache_parabolic
+  transform_variables!(u_transformed, u, equations_parabolic)
+
+  calc_gradient!(u_grad, u_transformed, t, mesh, equations_parabolic,
+                 boundary_conditions, dg, cache, cache_parabolic)
+
+  calc_viscous_fluxes!(viscous_flux, u_transformed, u_grad,
+                       mesh, equations_parabolic, dg, cache, cache_parabolic)
+
+  calc_divergence!(du, u_transformed, t, viscous_flux, mesh, equations_parabolic,
+                   boundary_conditions, dg, cache, cache_parabolic)
+
+  return nothing
+
+end
diff --git a/src/visualization/types.jl b/src/visualization/types.jl
index 4b4ab9cfeaf..392f120c461 100644
--- a/src/visualization/types.jl
+++ b/src/visualization/types.jl
@@ -406,6 +406,12 @@ struct ScalarData{T}
   data::T
 end
 
+"""
+    ScalarPlotData2D(u, semi::AbstractSemidiscretization; kwargs...)
+
+Returns an `PlotData2DTriangulated` object which is used to visualize a single scalar field.
+`u` should be an array whose entries correspond to values of the scalar field at nodal points.
+"""
 ScalarPlotData2D(u, semi::AbstractSemidiscretization; kwargs...) =
   ScalarPlotData2D(u, mesh_equations_solver_cache(semi)...; kwargs...)
 
diff --git a/test/runtests.jl b/test/runtests.jl
index e294d10746b..79a8e54e791 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -86,6 +86,10 @@ const TRIXI_NTHREADS   = clamp(Sys.CPU_THREADS, 2, 3)
     include("test_dgmulti_3d.jl")
   end
 
+  @time if TRIXI_TEST == "all" || TRIXI_TEST == "parabolic"
+    include("test_parabolic_2d.jl")
+  end
+
   @time if TRIXI_TEST == "all" || TRIXI_TEST == "misc_part1"
     include("test_unit.jl")
     include("test_visualization.jl")
diff --git a/test/test_dgmulti_2d.jl b/test/test_dgmulti_2d.jl
index d406be33c3e..90dd082c3ed 100644
--- a/test/test_dgmulti_2d.jl
+++ b/test/test_dgmulti_2d.jl
@@ -12,6 +12,7 @@ outdir = "out"
 isdir(outdir) && rm(outdir, recursive=true)
 
 @testset "DGMulti 2D" begin
+
   @trixi_testset "elixir_euler_weakform.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_weakform.jl"),
       cells_per_dimension = (4, 4),
diff --git a/test/test_parabolic_2d.jl b/test/test_parabolic_2d.jl
new file mode 100644
index 00000000000..8da8468c8c0
--- /dev/null
+++ b/test/test_parabolic_2d.jl
@@ -0,0 +1,83 @@
+module TestExamplesParabolic2D
+
+using Test
+using Trixi
+
+include("test_trixi.jl")
+
+EXAMPLES_DIR = joinpath(examples_dir(), "dgmulti_2d")
+
+# Start with a clean environment: remove Trixi output directory if it exists
+outdir = "out"
+isdir(outdir) && rm(outdir, recursive=true)
+
+@testset "SemidiscretizationHyperbolicParabolic" begin
+
+  @trixi_testset "DGMulti 2D rhs_parabolic!" begin
+
+    dg = DGMulti(polydeg = 2, element_type = Quad(), approximation_type = Polynomial(),
+                 surface_integral = SurfaceIntegralWeakForm(flux_central),
+                 volume_integral = VolumeIntegralWeakForm())
+    mesh = DGMultiMesh(dg, cells_per_dimension=(2, 2))
+
+    # test with polynomial initial condition x^2 * y
+    # test if we recover the exact second derivative
+    initial_condition = (x, t, equations) -> SVector(x[1]^2 * x[2])
+
+    equations = LinearScalarAdvectionEquation2D(1.0, 1.0)
+    equations_parabolic = LaplaceDiffusion2D(1.0, equations)
+
+    semi = SemidiscretizationHyperbolicParabolic(mesh, equations, equations_parabolic, initial_condition, dg)
+    @test_nowarn_debug show(stdout, semi)
+    @test_nowarn_debug show(stdout, MIME"text/plain"(), semi)
+
+    @test nvariables(semi)==nvariables(equations)
+    @test Base.ndims(semi)==Base.ndims(mesh)
+    @test Base.real(semi)==Base.real(dg)
+
+    ode = semidiscretize(semi, (0.0, 0.01))
+    u0 = similar(ode.u0)
+    Trixi.compute_coefficients!(u0, 0.0, semi)
+    @test u0 ≈ ode.u0
+
+    @unpack cache, cache_parabolic, equations_parabolic = semi
+    @unpack u_grad = cache_parabolic
+    for dim in eachindex(u_grad)
+      fill!(u_grad[dim], zero(eltype(u_grad[dim])))
+    end
+
+    t = 0.0
+    # pass in `boundary_condition_periodic` to fake "do-nothing"
+    # TODO: DGMulti. Make `boundary_condition_do_nothing` a callable struct like BoundaryConditionPeriodic
+    Trixi.calc_gradient!(u_grad, ode.u0, t, mesh, equations_parabolic,
+                         boundary_condition_periodic, dg, cache, cache_parabolic)
+    @unpack x, y = mesh.md
+    @test getindex.(u_grad[1], 1) ≈ 2 * x .* y
+    @test getindex.(u_grad[2], 1) ≈ x.^2
+
+    u_flux = similar.(u_grad)
+    Trixi.calc_viscous_fluxes!(u_flux, ode.u0, u_grad, mesh, equations_parabolic,
+                               dg, cache, cache_parabolic)
+    @test u_flux[1] ≈ u_grad[1]
+    @test u_flux[2] ≈ u_grad[2]
+
+    du = similar(ode.u0)
+    Trixi.calc_divergence!(du, ode.u0, t, u_flux, mesh, equations_parabolic,
+                           boundary_condition_periodic, dg, cache, cache_parabolic)
+    @test getindex.(du, 1) ≈ 2 * y
+  end
+
+  @trixi_testset "elixir_advection_diffusion.jl" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_diffusion.jl"),
+      cells_per_dimension = (4, 4), tspan=(0.0, 0.1),
+      l2 = [0.2485803335154642],
+      linf = [1.079606969242132]
+    )
+  end
+
+end
+
+# Clean up afterwards: delete Trixi output directory
+@test_nowarn isdir(outdir) && rm(outdir, recursive=true)
+
+end # module