diff --git a/examples/tree_2d_dgsem/elixir_navierstokes_taylor_green_vortex.jl b/examples/tree_2d_dgsem/elixir_navierstokes_taylor_green_vortex.jl
new file mode 100644
index 0000000000..c3cbc858f7
--- /dev/null
+++ b/examples/tree_2d_dgsem/elixir_navierstokes_taylor_green_vortex.jl
@@ -0,0 +1,78 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the compressible Navier-Stokes equations
+
+# TODO: parabolic; unify names of these accessor functions
+prandtl_number() = 0.72
+mu() = 6.25e-4 # equivalent to Re = 1600
+
+equations = CompressibleEulerEquations2D(1.4)
+equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, mu=mu(),
+                                                          Prandtl=prandtl_number())
+
+"""
+    initial_condition_taylor_green_vortex(x, t, equations::CompressibleEulerEquations2D)
+
+The classical viscous Taylor-Green vortex in 2D.
+This forms the basis behind the 3D case found for instance in
+  - Jonathan R. Bull and Antony Jameson
+  Simulation of the Compressible Taylor Green Vortex using High-Order Flux Reconstruction Schemes
+  [DOI: 10.2514/6.2014-3210](https://doi.org/10.2514/6.2014-3210)
+"""
+function initial_condition_taylor_green_vortex(x, t, equations::CompressibleEulerEquations2D)
+  A  = 1.0 # magnitude of speed
+  Ms = 0.1 # maximum Mach number
+
+  rho = 1.0
+  v1  =  A * sin(x[1]) * cos(x[2])
+  v2  = -A * cos(x[1]) * sin(x[2])
+  p   = (A / Ms)^2 * rho / equations.gamma # scaling to get Ms
+  p   = p + 1.0/4.0 * A^2 * rho * (cos(2*x[1]) + cos(2*x[2]))
+
+  return prim2cons(SVector(rho, v1, v2, p), equations)
+end
+initial_condition = initial_condition_taylor_green_vortex
+
+volume_flux = flux_ranocha
+solver = DGSEM(polydeg=3, surface_flux=flux_hllc,
+               volume_integral=VolumeIntegralFluxDifferencing(volume_flux))
+
+coordinates_min = (-1.0, -1.0) .* pi
+coordinates_max = ( 1.0,  1.0) .* pi
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level=4,
+                n_cells_max=100_000)
+
+
+semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic),
+                                             initial_condition, solver)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 20.0)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 100
+analysis_callback = AnalysisCallback(semi, interval=analysis_interval, save_analysis=true,
+                                     extra_analysis_integrals=(energy_kinetic,
+                                                               energy_internal))
+
+alive_callback = AliveCallback(analysis_interval=analysis_interval,)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback,
+                        alive_callback)
+
+###############################################################################
+# run the simulation
+
+time_int_tol = 1e-9
+sol = solve(ode, RDPK3SpFSAL49(); abstol=time_int_tol, reltol=time_int_tol,
+            ode_default_options()..., callback=callbacks)
+summary_callback() # print the timer summary
\ No newline at end of file
diff --git a/examples/tree_3d_dgsem/elixir_navierstokes_taylor_green_vortex.jl b/examples/tree_3d_dgsem/elixir_navierstokes_taylor_green_vortex.jl
index 9cb73a462b..5556831a59 100644
--- a/examples/tree_3d_dgsem/elixir_navierstokes_taylor_green_vortex.jl
+++ b/examples/tree_3d_dgsem/elixir_navierstokes_taylor_green_vortex.jl
@@ -16,7 +16,11 @@ equations_parabolic = CompressibleNavierStokesDiffusion3D(equations, mu=mu(),
 """
     initial_condition_taylor_green_vortex(x, t, equations::CompressibleEulerEquations3D)
 
-The classical inviscid Taylor-Green vortex.
+The classical viscous Taylor-Green vortex, as found for instance in
+
+- Jonathan R. Bull and Antony Jameson
+  Simulation of the Compressible Taylor Green Vortex using High-Order Flux Reconstruction Schemes
+  [DOI: 10.2514/6.2014-3210](https://doi.org/10.2514/6.2014-3210)
 """
 function initial_condition_taylor_green_vortex(x, t, equations::CompressibleEulerEquations3D)
   A  = 1.0 # magnitude of speed
diff --git a/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl b/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl
index 8f1e38c891..b12ecadb58 100644
--- a/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl
+++ b/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl
@@ -330,4 +330,19 @@ function rhs_parabolic!(du_ode, u_ode, semi::SemidiscretizationHyperbolicParabol
 
     return nothing
 end
+
+function _jacobian_ad_forward(semi::SemidiscretizationHyperbolicParabolic, t0, u0_ode,
+                              du_ode, config)
+    new_semi = remake(semi, uEltype = eltype(config))
+
+    du_ode_hyp = Vector{eltype(config)}(undef, length(du_ode))
+    J = ForwardDiff.jacobian(du_ode, u0_ode, config) do du_ode, u_ode
+        # Implementation of split ODE problem in OrdinaryDiffEq
+        rhs!(du_ode_hyp, u_ode, new_semi, t0)
+        rhs_parabolic!(du_ode, u_ode, new_semi, t0)
+        du_ode .+= du_ode_hyp
+    end
+
+    return J
+end
 end # @muladd
diff --git a/test/test_parabolic_2d.jl b/test/test_parabolic_2d.jl
index 57f296b55f..e3bb1ed9fb 100644
--- a/test/test_parabolic_2d.jl
+++ b/test/test_parabolic_2d.jl
@@ -188,6 +188,13 @@ isdir(outdir) && rm(outdir, recursive=true)
     )
   end
 
+  @trixi_testset "TreeMesh2D: elixir_navierstokes_taylor_green_vortex.jl" begin
+    @test_trixi_include(joinpath(examples_dir(), "tree_2d_dgsem", "elixir_navierstokes_taylor_green_vortex.jl"),
+      l2 = [0.0009279657228109691, 0.012454661988687185, 0.012454661988689886, 0.030487112728612178],
+      linf = [0.002435582543096171, 0.024824039368199546, 0.024824039368212758, 0.06731583711777489]
+    )
+  end
+
   @trixi_testset "P4estMesh2D: elixir_advection_diffusion_periodic.jl" begin
     @test_trixi_include(joinpath(examples_dir(), "p4est_2d_dgsem", "elixir_advection_diffusion_periodic.jl"),
       trees_per_dimension = (1, 1), initial_refinement_level = 2, tspan=(0.0, 0.5),
diff --git a/test/test_special_elixirs.jl b/test/test_special_elixirs.jl
index 23017059ea..c05dfbdfca 100644
--- a/test/test_special_elixirs.jl
+++ b/test/test_special_elixirs.jl
@@ -107,6 +107,15 @@ coverage = occursin("--code-coverage", cmd) && !occursin("--code-coverage=none",
       @test maximum(real, λ) < 10 * sqrt(eps(real(semi)))
     end
 
+    @timed_testset "Linear advection-diffusion" begin
+      trixi_include(@__MODULE__, joinpath(EXAMPLES_DIR, "tree_2d_dgsem", "elixir_advection_diffusion.jl"),
+                    tspan=(0.0, 0.0), initial_refinement_level=2)
+
+      J = jacobian_ad_forward(semi)
+      λ = eigvals(J)
+      @test maximum(real, λ) < 10 * sqrt(eps(real(semi)))
+    end
+
     @timed_testset "Compressible Euler equations" begin
       trixi_include(@__MODULE__, joinpath(EXAMPLES_DIR, "tree_2d_dgsem", "elixir_euler_density_wave.jl"),
                     tspan=(0.0, 0.0), initial_refinement_level=1)
@@ -165,6 +174,15 @@ coverage = occursin("--code-coverage", cmd) && !occursin("--code-coverage=none",
       end
     end
 
+    @timed_testset "Navier-Stokes" begin
+      trixi_include(@__MODULE__, joinpath(EXAMPLES_DIR, "tree_2d_dgsem", "elixir_navierstokes_taylor_green_vortex.jl"),
+                    tspan=(0.0, 0.0), initial_refinement_level=2)
+
+      J = jacobian_ad_forward(semi)
+      λ = eigvals(J)
+      @test maximum(real, λ) < 0.2
+    end
+
     @timed_testset "MHD" begin
       trixi_include(@__MODULE__, joinpath(EXAMPLES_DIR, "tree_2d_dgsem", "elixir_mhd_alfven_wave.jl"),
                     tspan=(0.0, 0.0), initial_refinement_level=0)