From d86a0f47bebe0945b3a70143015adb0fe95e6174 Mon Sep 17 00:00:00 2001
From: Benjamin Bolm <74359358+bennibolm@users.noreply.github.com>
Date: Tue, 26 Mar 2024 08:22:48 +0100
Subject: [PATCH 1/5] Use @batch reduction functionality for subcell bounds
 check (#1888)

* Use @batch reduction functionality for bounds check

* fmt

* Adapt compat bound for Polyester.jl

* Add comments
---
 Project.toml                                  |  2 +-
 docs/src/performance.md                       | 11 ----
 .../subcell_bounds_check_2d.jl                | 66 ++++++++-----------
 src/solvers/dgsem_tree/subcell_limiters_2d.jl | 11 +---
 4 files changed, 32 insertions(+), 58 deletions(-)

diff --git a/Project.toml b/Project.toml
index 27df49ed4fa..3db90c9928e 100644
--- a/Project.toml
+++ b/Project.toml
@@ -75,7 +75,7 @@ MuladdMacro = "0.2.2"
 Octavian = "0.3.21"
 OffsetArrays = "1.12"
 P4est = "0.4.9"
-Polyester = "0.7.5"
+Polyester = "0.7.10"
 PrecompileTools = "1.1"
 Preferences = "1.3"
 Printf = "1"
diff --git a/docs/src/performance.md b/docs/src/performance.md
index 40970e58c5c..9f81d3c3d8e 100644
--- a/docs/src/performance.md
+++ b/docs/src/performance.md
@@ -282,14 +282,3 @@ requires. It can thus be seen as a proxy for "energy used" and, as an extension,
     timing result, you need to set the analysis interval such that the
     `AnalysisCallback` is invoked at least once during the course of the simulation and
     discard the first PID value.
-
-## Performance issues with multi-threaded reductions
-[False sharing](https://en.wikipedia.org/wiki/False_sharing) is a known performance issue
-for systems with distributed caches. It also occurred for the implementation of a thread
-parallel bounds checking routine for the subcell IDP limiting
-in [PR #1736](https://github.com/trixi-framework/Trixi.jl/pull/1736).
-After some [testing and discussion](https://github.com/trixi-framework/Trixi.jl/pull/1736#discussion_r1423881895),
-it turned out that initializing a vector of length `n * Threads.nthreads()` and only using every
-n-th entry instead of a vector of length `Threads.nthreads()` fixes the problem.
-Since there are no processors with caches over 128B, we use `n = 128B / size(uEltype)`.
-Now, the bounds checking routine of the IDP limiting scales as hoped.
diff --git a/src/callbacks_stage/subcell_bounds_check_2d.jl b/src/callbacks_stage/subcell_bounds_check_2d.jl
index 19d73968c9a..3a56ea71f62 100644
--- a/src/callbacks_stage/subcell_bounds_check_2d.jl
+++ b/src/callbacks_stage/subcell_bounds_check_2d.jl
@@ -12,35 +12,37 @@
     (; variable_bounds) = limiter.cache.subcell_limiter_coefficients
     (; idp_bounds_delta_local, idp_bounds_delta_global) = limiter.cache
 
-    # Note: Accessing the threaded memory vector `idp_bounds_delta_local` with
-    # `deviation = idp_bounds_delta_local[key][Threads.threadid()]` causes critical performance
-    # issues due to False Sharing.
-    # Initializing a vector with n times the length and using every n-th entry fixes this
-    # problem and allows proper scaling:
-    # `deviation = idp_bounds_delta_local[key][n * Threads.threadid()]`
-    # Since there are no processors with caches over 128B, we use `n = 128B / size(uEltype)`
-    stride_size = div(128, sizeof(eltype(u))) # = n
+    # Note: In order to get the maximum deviation from the target bounds, this bounds check
+    # requires a reduction in every RK stage and for every enabled limiting option. To make
+    # this Thread-parallel we are using Polyester.jl's (at least v0.7.10) `@batch reduction`
+    # functionality.
+    # Although `@threaded` and `@batch` are currently used equivalently in Trixi.jl, we use
+    # `@batch` here to allow a possible redefinition of `@threaded` without creating errors here.
+    # See also https://github.com/trixi-framework/Trixi.jl/pull/1888#discussion_r1537785293.
 
     if local_minmax
         for v in limiter.local_minmax_variables_cons
             v_string = string(v)
             key_min = Symbol(v_string, "_min")
             key_max = Symbol(v_string, "_max")
-            deviation_min_threaded = idp_bounds_delta_local[key_min]
-            deviation_max_threaded = idp_bounds_delta_local[key_max]
-            @threaded for element in eachelement(solver, cache)
-                deviation_min = deviation_min_threaded[stride_size * Threads.threadid()]
-                deviation_max = deviation_max_threaded[stride_size * Threads.threadid()]
+            deviation_min = idp_bounds_delta_local[key_min]
+            deviation_max = idp_bounds_delta_local[key_max]
+            @batch reduction=((max, deviation_min), (max, deviation_max)) for element in eachelement(solver,
+                                                                                                     cache)
                 for j in eachnode(solver), i in eachnode(solver)
                     var = u[v, i, j, element]
+                    # Note: We always save the absolute deviations >= 0 and therefore use the
+                    # `max` operator for the lower and upper bound. The different directions of
+                    # upper and lower bound are considered in their calculations with a
+                    # different sign.
                     deviation_min = max(deviation_min,
                                         variable_bounds[key_min][i, j, element] - var)
                     deviation_max = max(deviation_max,
                                         var - variable_bounds[key_max][i, j, element])
                 end
-                deviation_min_threaded[stride_size * Threads.threadid()] = deviation_min
-                deviation_max_threaded[stride_size * Threads.threadid()] = deviation_max
             end
+            idp_bounds_delta_local[key_min] = deviation_min
+            idp_bounds_delta_local[key_max] = deviation_max
         end
     end
     if positivity
@@ -49,40 +51,35 @@
                 continue
             end
             key = Symbol(string(v), "_min")
-            deviation_threaded = idp_bounds_delta_local[key]
-            @threaded for element in eachelement(solver, cache)
-                deviation = deviation_threaded[stride_size * Threads.threadid()]
+            deviation = idp_bounds_delta_local[key]
+            @batch reduction=(max, deviation) for element in eachelement(solver, cache)
                 for j in eachnode(solver), i in eachnode(solver)
                     var = u[v, i, j, element]
                     deviation = max(deviation,
                                     variable_bounds[key][i, j, element] - var)
                 end
-                deviation_threaded[stride_size * Threads.threadid()] = deviation
             end
+            idp_bounds_delta_local[key] = deviation
         end
         for variable in limiter.positivity_variables_nonlinear
             key = Symbol(string(variable), "_min")
-            deviation_threaded = idp_bounds_delta_local[key]
-            @threaded for element in eachelement(solver, cache)
-                deviation = deviation_threaded[stride_size * Threads.threadid()]
+            deviation = idp_bounds_delta_local[key]
+            @batch reduction=(max, deviation) for element in eachelement(solver, cache)
                 for j in eachnode(solver), i in eachnode(solver)
                     var = variable(get_node_vars(u, equations, solver, i, j, element),
                                    equations)
                     deviation = max(deviation,
                                     variable_bounds[key][i, j, element] - var)
                 end
-                deviation_threaded[stride_size * Threads.threadid()] = deviation
             end
+            idp_bounds_delta_local[key] = deviation
         end
     end
 
     for (key, _) in idp_bounds_delta_local
-        # Calculate maximum deviations of all threads
-        idp_bounds_delta_local[key][stride_size] = maximum(idp_bounds_delta_local[key][stride_size * i]
-                                                           for i in 1:Threads.nthreads())
         # Update global maximum deviations
         idp_bounds_delta_global[key] = max(idp_bounds_delta_global[key],
-                                           idp_bounds_delta_local[key][stride_size])
+                                           idp_bounds_delta_local[key])
     end
 
     if save_errors
@@ -92,10 +89,8 @@
             if local_minmax
                 for v in limiter.local_minmax_variables_cons
                     v_string = string(v)
-                    print(f, ", ",
-                          idp_bounds_delta_local[Symbol(v_string, "_min")][stride_size],
-                          ", ",
-                          idp_bounds_delta_local[Symbol(v_string, "_max")][stride_size])
+                    print(f, ", ", idp_bounds_delta_local[Symbol(v_string, "_min")],
+                          ", ", idp_bounds_delta_local[Symbol(v_string, "_max")])
                 end
             end
             if positivity
@@ -103,21 +98,18 @@
                     if v in limiter.local_minmax_variables_cons
                         continue
                     end
-                    print(f, ", ",
-                          idp_bounds_delta_local[Symbol(string(v), "_min")][stride_size])
+                    print(f, ", ", idp_bounds_delta_local[Symbol(string(v), "_min")])
                 end
                 for variable in limiter.positivity_variables_nonlinear
                     print(f, ", ",
-                          idp_bounds_delta_local[Symbol(string(variable), "_min")][stride_size])
+                          idp_bounds_delta_local[Symbol(string(variable), "_min")])
                 end
             end
             println(f)
         end
         # Reset local maximum deviations
         for (key, _) in idp_bounds_delta_local
-            for i in 1:Threads.nthreads()
-                idp_bounds_delta_local[key][stride_size * i] = zero(eltype(idp_bounds_delta_local[key][stride_size]))
-            end
+            idp_bounds_delta_local[key] = zero(eltype(idp_bounds_delta_local[key]))
         end
     end
 
diff --git a/src/solvers/dgsem_tree/subcell_limiters_2d.jl b/src/solvers/dgsem_tree/subcell_limiters_2d.jl
index 3f7954c8958..9343cee4397 100644
--- a/src/solvers/dgsem_tree/subcell_limiters_2d.jl
+++ b/src/solvers/dgsem_tree/subcell_limiters_2d.jl
@@ -18,18 +18,11 @@ function create_cache(limiter::Type{SubcellLimiterIDP}, equations::AbstractEquat
 
     # Memory for bounds checking routine with `BoundsCheckCallback`.
     # Local variable contains the maximum deviation since the last export.
-    # Using a threaded vector to parallelize bounds check.
-    idp_bounds_delta_local = Dict{Symbol, Vector{real(basis)}}()
+    idp_bounds_delta_local = Dict{Symbol, real(basis)}()
     # Global variable contains the total maximum deviation.
     idp_bounds_delta_global = Dict{Symbol, real(basis)}()
-    # Note: False sharing causes critical performance issues on multiple threads when using a vector
-    # of length `Threads.nthreads()`. Initializing a vector of length `n * Threads.nthreads()`
-    # and then only using every n-th entry, fixes the problem and allows proper scaling.
-    # Since there are no processors with caches over 128B, we use `n = 128B / size(uEltype)`
-    stride_size = div(128, sizeof(eltype(basis.nodes))) # = n
     for key in bound_keys
-        idp_bounds_delta_local[key] = [zero(real(basis))
-                                       for _ in 1:(stride_size * Threads.nthreads())]
+        idp_bounds_delta_local[key] = zero(real(basis))
         idp_bounds_delta_global[key] = zero(real(basis))
     end
 

From cd6fdaaa1443376ad97bc408cafc27d82a647175 Mon Sep 17 00:00:00 2001
From: Benedict <135045760+bgeihe@users.noreply.github.com>
Date: Tue, 26 Mar 2024 08:23:15 +0100
Subject: [PATCH 2/5] abort initial AMR loop after 10 iterations (#1890)

Co-authored-by: Hendrik Ranocha <ranocha@users.noreply.github.com>
---
 src/callbacks_step/amr.jl | 7 +++++++
 1 file changed, 7 insertions(+)

diff --git a/src/callbacks_step/amr.jl b/src/callbacks_step/amr.jl
index 6f57d6647fc..1ab65a3553e 100644
--- a/src/callbacks_step/amr.jl
+++ b/src/callbacks_step/amr.jl
@@ -138,11 +138,18 @@ function initialize!(cb::DiscreteCallback{Condition, Affect!}, u, t,
         # iterate until mesh does not change anymore
         has_changed = amr_callback(integrator,
                                    only_refine = amr_callback.adapt_initial_condition_only_refine)
+        iterations = 1
         while has_changed
             compute_coefficients!(integrator.u, t, semi)
             u_modified!(integrator, true)
             has_changed = amr_callback(integrator,
                                        only_refine = amr_callback.adapt_initial_condition_only_refine)
+            iterations = iterations + 1
+            if iterations > 10
+                @warn "AMR for initial condition did not settle within 10 iterations!\n" *
+                      "Consider adjusting thresholds or setting `adapt_initial_condition_only_refine`."
+                break
+            end
         end
     end
 

From 27026de6e2be4f0f6655efb68c082515f1383e84 Mon Sep 17 00:00:00 2001
From: Hendrik Ranocha <mail@ranocha.de>
Date: Tue, 26 Mar 2024 10:07:11 +0100
Subject: [PATCH 3/5] set version to v0.7.5

---
 Project.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index 3db90c9928e..9ca3086306f 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "Trixi"
 uuid = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb"
 authors = ["Michael Schlottke-Lakemper <m.schlottke-lakemper@hlrs.de>", "Gregor Gassner <ggassner@uni-koeln.de>", "Hendrik Ranocha <mail@ranocha.de>", "Andrew R. Winters <andrew.ross.winters@liu.se>", "Jesse Chan <jesse.chan@rice.edu>"]
-version = "0.7.5-pre"
+version = "0.7.5"
 
 [deps]
 CodeTracking = "da1fd8a2-8d9e-5ec2-8556-3022fb5608a2"

From 909abb4473c95042a87e93e46d443dc381c2b523 Mon Sep 17 00:00:00 2001
From: Hendrik Ranocha <mail@ranocha.de>
Date: Tue, 26 Mar 2024 10:07:23 +0100
Subject: [PATCH 4/5] set development version to v0.7.6-pre

---
 Project.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index 9ca3086306f..6ff7f29686d 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "Trixi"
 uuid = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb"
 authors = ["Michael Schlottke-Lakemper <m.schlottke-lakemper@hlrs.de>", "Gregor Gassner <ggassner@uni-koeln.de>", "Hendrik Ranocha <mail@ranocha.de>", "Andrew R. Winters <andrew.ross.winters@liu.se>", "Jesse Chan <jesse.chan@rice.edu>"]
-version = "0.7.5"
+version = "0.7.6-pre"
 
 [deps]
 CodeTracking = "da1fd8a2-8d9e-5ec2-8556-3022fb5608a2"

From 73da92c446c6f60fc09ad1cf5e579bbf5eccb34d Mon Sep 17 00:00:00 2001
From: Daniel Doehring <daniel.doehring@rwth-aachen.de>
Date: Wed, 27 Mar 2024 15:55:37 +0100
Subject: [PATCH 5/5] Sutherlands Law for temperature dependent viscosity
 (#1808)

* sample 2D

* bug fixes

* typo

* sutherlands 1d

* sutherlands 3d

* fmt

* main merge

* test var mu

* variable mu

* fmt

* example using sutherland

* comment

* example

* reset toml

* Shorten

* Unit test

* fmt

* Update examples/tree_2d_dgsem/elixir_navierstokes_taylor_green_vortex_sutherland.jl

* remove todo

* Apply suggestions from code review

Co-authored-by: Andrew Winters <andrew.ross.winters@liu.se>

* test vals

* non-functionalized mu

* comments & dosctrings

* docstrings

* fix tests

* avoid if

* docstring and comments

* Update src/equations/compressible_navier_stokes.jl

Co-authored-by: Michael Schlottke-Lakemper <michael@sloede.com>

* fmt

---------

Co-authored-by: Andrew Winters <andrew.ross.winters@liu.se>
Co-authored-by: Michael Schlottke-Lakemper <michael@sloede.com>
Co-authored-by: Jesse Chan <1156048+jlchan@users.noreply.github.com>
---
 .../elixir_navierstokes_lid_driven_cavity.jl  |  5 +-
 ...elixir_navierstokes_taylor_green_vortex.jl |  5 +-
 .../elixir_navierstokes_lid_driven_cavity.jl  |  5 +-
 ...ixir_navierstokes_lid_driven_cavity_amr.jl |  5 +-
 .../elixir_navierstokes_blast_wave_amr.jl     |  5 +-
 ...elixir_navierstokes_taylor_green_vortex.jl |  5 +-
 ...ir_navierstokes_taylor_green_vortex_amr.jl |  5 +-
 ...lixir_navierstokes_convergence_periodic.jl |  1 -
 .../elixir_navierstokes_lid_driven_cavity.jl  |  5 +-
 .../elixir_navierstokes_shearlayer_amr.jl     |  5 +-
 ...elixir_navierstokes_taylor_green_vortex.jl |  5 +-
 ...erstokes_taylor_green_vortex_sutherland.jl | 94 +++++++++++++++++++
 ...elixir_navierstokes_taylor_green_vortex.jl |  5 +-
 src/equations/compressible_navier_stokes.jl   | 14 +++
 .../compressible_navier_stokes_1d.jl          | 26 +++--
 .../compressible_navier_stokes_2d.jl          | 26 +++--
 .../compressible_navier_stokes_3d.jl          | 26 +++--
 test/test_parabolic_2d.jl                     | 34 +++++--
 test/test_unit.jl                             | 39 ++++++++
 19 files changed, 243 insertions(+), 72 deletions(-)
 create mode 100644 examples/tree_2d_dgsem/elixir_navierstokes_taylor_green_vortex_sutherland.jl

diff --git a/examples/dgmulti_2d/elixir_navierstokes_lid_driven_cavity.jl b/examples/dgmulti_2d/elixir_navierstokes_lid_driven_cavity.jl
index 7c55cbf0ccf..a612dd0e0d5 100644
--- a/examples/dgmulti_2d/elixir_navierstokes_lid_driven_cavity.jl
+++ b/examples/dgmulti_2d/elixir_navierstokes_lid_driven_cavity.jl
@@ -4,12 +4,11 @@ using Trixi
 ###############################################################################
 # semidiscretization of the ideal compressible Navier-Stokes equations
 
-# TODO: parabolic; unify names of these accessor functions
 prandtl_number() = 0.72
-mu() = 0.001
+mu = 0.001
 
 equations = CompressibleEulerEquations2D(1.4)
-equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, mu = mu(),
+equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, mu = mu,
                                                           Prandtl = prandtl_number())
 
 # Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux
diff --git a/examples/dgmulti_3d/elixir_navierstokes_taylor_green_vortex.jl b/examples/dgmulti_3d/elixir_navierstokes_taylor_green_vortex.jl
index dedd8267a3b..9ae90ac47b6 100644
--- a/examples/dgmulti_3d/elixir_navierstokes_taylor_green_vortex.jl
+++ b/examples/dgmulti_3d/elixir_navierstokes_taylor_green_vortex.jl
@@ -5,12 +5,11 @@ 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
+mu = 6.25e-4 # equivalent to Re = 1600
 
 equations = CompressibleEulerEquations3D(1.4)
-equations_parabolic = CompressibleNavierStokesDiffusion3D(equations, mu = mu(),
+equations_parabolic = CompressibleNavierStokesDiffusion3D(equations, mu = mu,
                                                           Prandtl = prandtl_number())
 
 """
diff --git a/examples/p4est_2d_dgsem/elixir_navierstokes_lid_driven_cavity.jl b/examples/p4est_2d_dgsem/elixir_navierstokes_lid_driven_cavity.jl
index bc28ae6ffb3..728736fe49e 100644
--- a/examples/p4est_2d_dgsem/elixir_navierstokes_lid_driven_cavity.jl
+++ b/examples/p4est_2d_dgsem/elixir_navierstokes_lid_driven_cavity.jl
@@ -4,12 +4,11 @@ using Trixi
 ###############################################################################
 # semidiscretization of the ideal compressible Navier-Stokes equations
 
-# TODO: parabolic; unify names of these accessor functions
 prandtl_number() = 0.72
-mu() = 0.001
+mu = 0.001
 
 equations = CompressibleEulerEquations2D(1.4)
-equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, mu = mu(),
+equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, mu = mu,
                                                           Prandtl = prandtl_number())
 
 # Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux
diff --git a/examples/p4est_2d_dgsem/elixir_navierstokes_lid_driven_cavity_amr.jl b/examples/p4est_2d_dgsem/elixir_navierstokes_lid_driven_cavity_amr.jl
index 898366969a4..d9eaf4fdb72 100644
--- a/examples/p4est_2d_dgsem/elixir_navierstokes_lid_driven_cavity_amr.jl
+++ b/examples/p4est_2d_dgsem/elixir_navierstokes_lid_driven_cavity_amr.jl
@@ -4,12 +4,11 @@ using Trixi
 ###############################################################################
 # semidiscretization of the ideal compressible Navier-Stokes equations
 
-# TODO: parabolic; unify names of these accessor functions
 prandtl_number() = 0.72
-mu() = 0.001
+mu = 0.001
 
 equations = CompressibleEulerEquations2D(1.4)
-equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, mu = mu(),
+equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, mu = mu,
                                                           Prandtl = prandtl_number())
 
 # Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux
diff --git a/examples/p4est_3d_dgsem/elixir_navierstokes_blast_wave_amr.jl b/examples/p4est_3d_dgsem/elixir_navierstokes_blast_wave_amr.jl
index 5df89fbcdf2..d556d0ab70d 100644
--- a/examples/p4est_3d_dgsem/elixir_navierstokes_blast_wave_amr.jl
+++ b/examples/p4est_3d_dgsem/elixir_navierstokes_blast_wave_amr.jl
@@ -5,12 +5,11 @@ 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
+mu = 6.25e-4 # equivalent to Re = 1600
 
 equations = CompressibleEulerEquations3D(1.4)
-equations_parabolic = CompressibleNavierStokesDiffusion3D(equations, mu = mu(),
+equations_parabolic = CompressibleNavierStokesDiffusion3D(equations, mu = mu,
                                                           Prandtl = prandtl_number())
 
 function initial_condition_3d_blast_wave(x, t, equations::CompressibleEulerEquations3D)
diff --git a/examples/p4est_3d_dgsem/elixir_navierstokes_taylor_green_vortex.jl b/examples/p4est_3d_dgsem/elixir_navierstokes_taylor_green_vortex.jl
index d785013f5a9..9c90e4d3218 100644
--- a/examples/p4est_3d_dgsem/elixir_navierstokes_taylor_green_vortex.jl
+++ b/examples/p4est_3d_dgsem/elixir_navierstokes_taylor_green_vortex.jl
@@ -5,12 +5,11 @@ 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
+mu = 6.25e-4 # equivalent to Re = 1600
 
 equations = CompressibleEulerEquations3D(1.4)
-equations_parabolic = CompressibleNavierStokesDiffusion3D(equations, mu = mu(),
+equations_parabolic = CompressibleNavierStokesDiffusion3D(equations, mu = mu,
                                                           Prandtl = prandtl_number())
 
 """
diff --git a/examples/p4est_3d_dgsem/elixir_navierstokes_taylor_green_vortex_amr.jl b/examples/p4est_3d_dgsem/elixir_navierstokes_taylor_green_vortex_amr.jl
index c15227a1c29..2741f0df174 100644
--- a/examples/p4est_3d_dgsem/elixir_navierstokes_taylor_green_vortex_amr.jl
+++ b/examples/p4est_3d_dgsem/elixir_navierstokes_taylor_green_vortex_amr.jl
@@ -5,12 +5,11 @@ 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
+mu = 6.25e-4 # equivalent to Re = 1600
 
 equations = CompressibleEulerEquations3D(1.4)
-equations_parabolic = CompressibleNavierStokesDiffusion3D(equations, mu = mu(),
+equations_parabolic = CompressibleNavierStokesDiffusion3D(equations, mu = mu,
                                                           Prandtl = prandtl_number())
 
 """
diff --git a/examples/tree_1d_dgsem/elixir_navierstokes_convergence_periodic.jl b/examples/tree_1d_dgsem/elixir_navierstokes_convergence_periodic.jl
index 33ad0d5271c..eab0840f385 100644
--- a/examples/tree_1d_dgsem/elixir_navierstokes_convergence_periodic.jl
+++ b/examples/tree_1d_dgsem/elixir_navierstokes_convergence_periodic.jl
@@ -5,7 +5,6 @@ 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
 
diff --git a/examples/tree_2d_dgsem/elixir_navierstokes_lid_driven_cavity.jl b/examples/tree_2d_dgsem/elixir_navierstokes_lid_driven_cavity.jl
index 70d76fc9075..b8e20e27f66 100644
--- a/examples/tree_2d_dgsem/elixir_navierstokes_lid_driven_cavity.jl
+++ b/examples/tree_2d_dgsem/elixir_navierstokes_lid_driven_cavity.jl
@@ -4,12 +4,11 @@ using Trixi
 ###############################################################################
 # semidiscretization of the ideal compressible Navier-Stokes equations
 
-# TODO: parabolic; unify names of these accessor functions
 prandtl_number() = 0.72
-mu() = 0.001
+mu = 0.001
 
 equations = CompressibleEulerEquations2D(1.4)
-equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, mu = mu(),
+equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, mu = mu,
                                                           Prandtl = prandtl_number())
 
 # Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux
diff --git a/examples/tree_2d_dgsem/elixir_navierstokes_shearlayer_amr.jl b/examples/tree_2d_dgsem/elixir_navierstokes_shearlayer_amr.jl
index 4d92ea261e9..a7492bafb47 100644
--- a/examples/tree_2d_dgsem/elixir_navierstokes_shearlayer_amr.jl
+++ b/examples/tree_2d_dgsem/elixir_navierstokes_shearlayer_amr.jl
@@ -5,12 +5,11 @@ using Trixi
 ###############################################################################
 # semidiscretization of the compressible Navier-Stokes equations
 
-# TODO: parabolic; unify names of these accessor functions
 prandtl_number() = 0.72
-mu() = 1.0 / 3.0 * 10^(-5) # equivalent to Re = 30,000
+mu = 1.0 / 3.0 * 10^(-4) # equivalent to Re = 30,000
 
 equations = CompressibleEulerEquations2D(1.4)
-equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, mu = mu(),
+equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, mu = mu,
                                                           Prandtl = prandtl_number())
 
 """
diff --git a/examples/tree_2d_dgsem/elixir_navierstokes_taylor_green_vortex.jl b/examples/tree_2d_dgsem/elixir_navierstokes_taylor_green_vortex.jl
index a3e38bf6d92..c6e5f0bc40a 100644
--- a/examples/tree_2d_dgsem/elixir_navierstokes_taylor_green_vortex.jl
+++ b/examples/tree_2d_dgsem/elixir_navierstokes_taylor_green_vortex.jl
@@ -5,12 +5,11 @@ 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
+mu = 6.25e-4 # equivalent to Re = 1600
 
 equations = CompressibleEulerEquations2D(1.4)
-equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, mu = mu(),
+equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, mu = mu,
                                                           Prandtl = prandtl_number())
 
 """
diff --git a/examples/tree_2d_dgsem/elixir_navierstokes_taylor_green_vortex_sutherland.jl b/examples/tree_2d_dgsem/elixir_navierstokes_taylor_green_vortex_sutherland.jl
new file mode 100644
index 00000000000..9598ae184cd
--- /dev/null
+++ b/examples/tree_2d_dgsem/elixir_navierstokes_taylor_green_vortex_sutherland.jl
@@ -0,0 +1,94 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the compressible Navier-Stokes equations
+
+prandtl_number() = 0.72
+
+# Use Sutherland's law for a temperature-dependent viscosity.
+# For details, see e.g. 
+# Frank M. White: Viscous Fluid Flow, 2nd Edition.
+# 1991, McGraw-Hill, ISBN, 0-07-069712-4
+# Pages 28 and 29.
+@inline function mu(u, equations)
+    T_ref = 291.15
+
+    R_specific_air = 287.052874
+    T = R_specific_air * Trixi.temperature(u, equations)
+
+    C_air = 120.0
+    mu_ref_air = 1.827e-5
+
+    return mu_ref_air * (T_ref + C_air) / (T + C_air) * (T / T_ref)^1.5
+end
+
+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
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 65bd9aa133d..3e54c791ec6 100644
--- a/examples/tree_3d_dgsem/elixir_navierstokes_taylor_green_vortex.jl
+++ b/examples/tree_3d_dgsem/elixir_navierstokes_taylor_green_vortex.jl
@@ -5,12 +5,11 @@ 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
+mu = 6.25e-4 # equivalent to Re = 1600
 
 equations = CompressibleEulerEquations3D(1.4)
-equations_parabolic = CompressibleNavierStokesDiffusion3D(equations, mu = mu(),
+equations_parabolic = CompressibleNavierStokesDiffusion3D(equations, mu = mu,
                                                           Prandtl = prandtl_number())
 
 """
diff --git a/src/equations/compressible_navier_stokes.jl b/src/equations/compressible_navier_stokes.jl
index 3059771197c..c530625f226 100644
--- a/src/equations/compressible_navier_stokes.jl
+++ b/src/equations/compressible_navier_stokes.jl
@@ -62,3 +62,17 @@ Under `GradientVariablesEntropy`, the Navier-Stokes discretization is provably e
 """
 struct GradientVariablesPrimitive end
 struct GradientVariablesEntropy end
+
+"""
+    dynamic_viscosity(u, equations)
+
+Wrapper for the dynamic viscosity that calls
+`dynamic_viscosity(u, equations.mu, equations)`, which dispatches on the type of 
+`equations.mu`. 
+For constant `equations.mu`, i.e., `equations.mu` is of `Real`-type it is returned directly.
+In all other cases, `equations.mu` is assumed to be a function with arguments
+`u` and `equations` and is called with these arguments.
+"""
+dynamic_viscosity(u, equations) = dynamic_viscosity(u, equations.mu, equations)
+dynamic_viscosity(u, mu::Real, equations) = mu
+dynamic_viscosity(u, mu::T, equations) where {T} = mu(u, equations)
diff --git a/src/equations/compressible_navier_stokes_1d.jl b/src/equations/compressible_navier_stokes_1d.jl
index d2c46ecc7d8..3dbdf11c56b 100644
--- a/src/equations/compressible_navier_stokes_1d.jl
+++ b/src/equations/compressible_navier_stokes_1d.jl
@@ -21,6 +21,9 @@ the [`CompressibleEulerEquations1D`](@ref).
 
 Fluid properties such as the dynamic viscosity ``\mu`` can be provided in any consistent unit system, e.g.,
 [``\mu``] = kg m⁻¹ s⁻¹.
+The viscosity ``\mu`` may be a constant or a function of the current state, e.g., 
+depending on temperature (Sutherland's law): ``\mu = \mu(T)``.
+In the latter case, the function `mu` needs to have the signature `mu(u, equations)`.
 
 The particular form of the compressible Navier-Stokes implemented is
 ```math
@@ -80,7 +83,7 @@ where
 w_2 = \frac{\rho v1}{p},\, w_3 = -\frac{\rho}{p}
 ```
 """
-struct CompressibleNavierStokesDiffusion1D{GradientVariables, RealT <: Real,
+struct CompressibleNavierStokesDiffusion1D{GradientVariables, RealT <: Real, Mu,
                                            E <: AbstractCompressibleEulerEquations{1}} <:
        AbstractCompressibleNavierStokesDiffusion{1, 3, GradientVariables}
     # TODO: parabolic
@@ -89,7 +92,7 @@ struct CompressibleNavierStokesDiffusion1D{GradientVariables, RealT <: Real,
     gamma::RealT               # ratio of specific heats
     inv_gamma_minus_one::RealT # = inv(gamma - 1); can be used to write slow divisions as fast multiplications
 
-    mu::RealT                  # viscosity
+    mu::Mu                     # viscosity
     Pr::RealT                  # Prandtl number
     kappa::RealT               # thermal diffusivity for Fick's law
 
@@ -103,16 +106,17 @@ function CompressibleNavierStokesDiffusion1D(equations::CompressibleEulerEquatio
                                              gradient_variables = GradientVariablesPrimitive())
     gamma = equations.gamma
     inv_gamma_minus_one = equations.inv_gamma_minus_one
-    μ, Pr = promote(mu, Prandtl)
 
     # Under the assumption of constant Prandtl number the thermal conductivity
-    # constant is kappa = gamma μ / ((gamma-1) Pr).
+    # constant is kappa = gamma μ / ((gamma-1) Prandtl).
     # Important note! Factor of μ is accounted for later in `flux`.
-    kappa = gamma * inv_gamma_minus_one / Pr
+    # This avoids recomputation of kappa for non-constant μ.
+    kappa = gamma * inv_gamma_minus_one / Prandtl
 
     CompressibleNavierStokesDiffusion1D{typeof(gradient_variables), typeof(gamma),
+                                        typeof(mu),
                                         typeof(equations)}(gamma, inv_gamma_minus_one,
-                                                           μ, Pr, kappa,
+                                                           mu, Prandtl, kappa,
                                                            equations,
                                                            gradient_variables)
 end
@@ -159,10 +163,12 @@ function flux(u, gradients, orientation::Integer,
     # in the implementation
     q1 = equations.kappa * dTdx
 
-    # Constant dynamic viscosity is copied to a variable for readability.
-    # Offers flexibility for dynamic viscosity via Sutherland's law where it depends
-    # on temperature and reference values, Ts and Tref such that mu(T)
-    mu = equations.mu
+    # In the simplest cases, the user passed in `mu` or `mu()` 
+    # (which returns just a constant) but
+    # more complex functions like Sutherland's law are possible.
+    # `dynamic_viscosity` is a helper function that handles both cases
+    # by dispatching on the type of `equations.mu`.
+    mu = dynamic_viscosity(u, equations)
 
     # viscous flux components in the x-direction
     f1 = zero(rho)
diff --git a/src/equations/compressible_navier_stokes_2d.jl b/src/equations/compressible_navier_stokes_2d.jl
index 5df7c01ca5c..3256343703a 100644
--- a/src/equations/compressible_navier_stokes_2d.jl
+++ b/src/equations/compressible_navier_stokes_2d.jl
@@ -21,6 +21,9 @@ the [`CompressibleEulerEquations2D`](@ref).
 
 Fluid properties such as the dynamic viscosity ``\mu`` can be provided in any consistent unit system, e.g.,
 [``\mu``] = kg m⁻¹ s⁻¹.
+The viscosity ``\mu`` may be a constant or a function of the current state, e.g., 
+depending on temperature (Sutherland's law): ``\mu = \mu(T)``.
+In the latter case, the function `mu` needs to have the signature `mu(u, equations)`.
 
 The particular form of the compressible Navier-Stokes implemented is
 ```math
@@ -80,7 +83,7 @@ where
 w_2 = \frac{\rho v_1}{p},\, w_3 = \frac{\rho v_2}{p},\, w_4 = -\frac{\rho}{p}
 ```
 """
-struct CompressibleNavierStokesDiffusion2D{GradientVariables, RealT <: Real,
+struct CompressibleNavierStokesDiffusion2D{GradientVariables, RealT <: Real, Mu,
                                            E <: AbstractCompressibleEulerEquations{2}} <:
        AbstractCompressibleNavierStokesDiffusion{2, 4, GradientVariables}
     # TODO: parabolic
@@ -89,7 +92,7 @@ struct CompressibleNavierStokesDiffusion2D{GradientVariables, RealT <: Real,
     gamma::RealT               # ratio of specific heats
     inv_gamma_minus_one::RealT # = inv(gamma - 1); can be used to write slow divisions as fast multiplications
 
-    mu::RealT                  # viscosity
+    mu::Mu                     # viscosity
     Pr::RealT                  # Prandtl number
     kappa::RealT               # thermal diffusivity for Fick's law
 
@@ -103,16 +106,17 @@ function CompressibleNavierStokesDiffusion2D(equations::CompressibleEulerEquatio
                                              gradient_variables = GradientVariablesPrimitive())
     gamma = equations.gamma
     inv_gamma_minus_one = equations.inv_gamma_minus_one
-    μ, Pr = promote(mu, Prandtl)
 
     # Under the assumption of constant Prandtl number the thermal conductivity
-    # constant is kappa = gamma μ / ((gamma-1) Pr).
+    # constant is kappa = gamma μ / ((gamma-1) Prandtl).
     # Important note! Factor of μ is accounted for later in `flux`.
-    kappa = gamma * inv_gamma_minus_one / Pr
+    # This avoids recomputation of kappa for non-constant μ.
+    kappa = gamma * inv_gamma_minus_one / Prandtl
 
     CompressibleNavierStokesDiffusion2D{typeof(gradient_variables), typeof(gamma),
+                                        typeof(mu),
                                         typeof(equations)}(gamma, inv_gamma_minus_one,
-                                                           μ, Pr, kappa,
+                                                           mu, Prandtl, kappa,
                                                            equations,
                                                            gradient_variables)
 end
@@ -168,10 +172,12 @@ function flux(u, gradients, orientation::Integer,
     q1 = equations.kappa * dTdx
     q2 = equations.kappa * dTdy
 
-    # Constant dynamic viscosity is copied to a variable for readability.
-    # Offers flexibility for dynamic viscosity via Sutherland's law where it depends
-    # on temperature and reference values, Ts and Tref such that mu(T)
-    mu = equations.mu
+    # In the simplest cases, the user passed in `mu` or `mu()` 
+    # (which returns just a constant) but
+    # more complex functions like Sutherland's law are possible.
+    # `dynamic_viscosity` is a helper function that handles both cases
+    # by dispatching on the type of `equations.mu`.
+    mu = dynamic_viscosity(u, equations)
 
     if orientation == 1
         # viscous flux components in the x-direction
diff --git a/src/equations/compressible_navier_stokes_3d.jl b/src/equations/compressible_navier_stokes_3d.jl
index e5567ae5789..9833122eb32 100644
--- a/src/equations/compressible_navier_stokes_3d.jl
+++ b/src/equations/compressible_navier_stokes_3d.jl
@@ -21,6 +21,9 @@ the [`CompressibleEulerEquations3D`](@ref).
 
 Fluid properties such as the dynamic viscosity ``\mu`` can be provided in any consistent unit system, e.g.,
 [``\mu``] = kg m⁻¹ s⁻¹.
+The viscosity ``\mu`` may be a constant or a function of the current state, e.g., 
+depending on temperature (Sutherland's law): ``\mu = \mu(T)``.
+In the latter case, the function `mu` needs to have the signature `mu(u, equations)`.
 
 The particular form of the compressible Navier-Stokes implemented is
 ```math
@@ -80,7 +83,7 @@ where
 w_2 = \frac{\rho v_1}{p},\, w_3 = \frac{\rho v_2}{p},\, w_4 = \frac{\rho v_3}{p},\, w_5 = -\frac{\rho}{p}
 ```
 """
-struct CompressibleNavierStokesDiffusion3D{GradientVariables, RealT <: Real,
+struct CompressibleNavierStokesDiffusion3D{GradientVariables, RealT <: Real, Mu,
                                            E <: AbstractCompressibleEulerEquations{3}} <:
        AbstractCompressibleNavierStokesDiffusion{3, 5, GradientVariables}
     # TODO: parabolic
@@ -89,7 +92,7 @@ struct CompressibleNavierStokesDiffusion3D{GradientVariables, RealT <: Real,
     gamma::RealT               # ratio of specific heats
     inv_gamma_minus_one::RealT # = inv(gamma - 1); can be used to write slow divisions as fast multiplications
 
-    mu::RealT                  # viscosity
+    mu::Mu                     # viscosity
     Pr::RealT                  # Prandtl number
     kappa::RealT               # thermal diffusivity for Fick's law
 
@@ -103,16 +106,17 @@ function CompressibleNavierStokesDiffusion3D(equations::CompressibleEulerEquatio
                                              gradient_variables = GradientVariablesPrimitive())
     gamma = equations.gamma
     inv_gamma_minus_one = equations.inv_gamma_minus_one
-    μ, Pr = promote(mu, Prandtl)
 
     # Under the assumption of constant Prandtl number the thermal conductivity
-    # constant is kappa = gamma μ / ((gamma-1) Pr).
+    # constant is kappa = gamma μ / ((gamma-1) Prandtl).
     # Important note! Factor of μ is accounted for later in `flux`.
-    kappa = gamma * inv_gamma_minus_one / Pr
+    # This avoids recomputation of kappa for non-constant μ.
+    kappa = gamma * inv_gamma_minus_one / Prandtl
 
     CompressibleNavierStokesDiffusion3D{typeof(gradient_variables), typeof(gamma),
+                                        typeof(mu),
                                         typeof(equations)}(gamma, inv_gamma_minus_one,
-                                                           μ, Pr, kappa,
+                                                           mu, Prandtl, kappa,
                                                            equations,
                                                            gradient_variables)
 end
@@ -181,10 +185,12 @@ function flux(u, gradients, orientation::Integer,
     q2 = equations.kappa * dTdy
     q3 = equations.kappa * dTdz
 
-    # Constant dynamic viscosity is copied to a variable for readability.
-    # Offers flexibility for dynamic viscosity via Sutherland's law where it depends
-    # on temperature and reference values, Ts and Tref such that mu(T)
-    mu = equations.mu
+    # In the simplest cases, the user passed in `mu` or `mu()` 
+    # (which returns just a constant) but
+    # more complex functions like Sutherland's law are possible.
+    # `dynamic_viscosity` is a helper function that handles both cases
+    # by dispatching on the type of `equations.mu`.
+    mu = dynamic_viscosity(u, equations)
 
     if orientation == 1
         # viscous flux components in the x-direction
diff --git a/test/test_parabolic_2d.jl b/test/test_parabolic_2d.jl
index 9f1382caa62..d47c34f9e75 100644
--- a/test/test_parabolic_2d.jl
+++ b/test/test_parabolic_2d.jl
@@ -476,20 +476,38 @@ end
     @test_trixi_include(joinpath(examples_dir(), "tree_2d_dgsem",
                                  "elixir_navierstokes_shearlayer_amr.jl"),
                         l2=[
-                            0.00526017743452336,
-                            0.4130430692895672,
-                            0.4310996183791349,
-                            1.1544344171604635,
+                            0.005155557460409018,
+                            0.4048446934219344,
+                            0.43040068852937047,
+                            1.1255130552079322,
                         ],
                         linf=[
-                            0.03492185879198495,
-                            1.392635891671335,
-                            1.357551616406459,
-                            8.713760873018146,
+                            0.03287305649809613,
+                            1.1656793717431393,
+                            1.3917196016246969,
+                            8.146587380114653,
                         ],
                         tspan=(0.0, 0.7))
 end
 
+@trixi_testset "TreeMesh2D: elixir_navierstokes_taylor_green_vortex_sutherland.jl" begin
+    @test_trixi_include(joinpath(examples_dir(), "tree_2d_dgsem",
+                                 "elixir_navierstokes_taylor_green_vortex_sutherland.jl"),
+                        l2=[
+                            0.001452856280034929,
+                            0.0007538775539989481,
+                            0.0007538775539988681,
+                            0.011035506549989587,
+                        ],
+                        linf=[
+                            0.003291912841311362,
+                            0.002986462478096974,
+                            0.0029864624780958637,
+                            0.0231954665514138,
+                        ],
+                        tspan=(0.0, 1.0))
+end
+
 @trixi_testset "P4estMesh2D: elixir_advection_diffusion_periodic.jl" begin
     @test_trixi_include(joinpath(examples_dir(), "p4est_2d_dgsem",
                                  "elixir_advection_diffusion_periodic.jl"),
diff --git a/test/test_unit.jl b/test/test_unit.jl
index 03a78f6918a..79950f58d59 100644
--- a/test/test_unit.jl
+++ b/test/test_unit.jl
@@ -1576,6 +1576,45 @@ end
 
     @test mesh.boundary_faces[:entire_boundary] == [1, 2]
 end
+
+@testset "Sutherlands Law" begin
+    function mu(u, equations)
+        T_ref = 291.15
+
+        R_specific_air = 287.052874
+        T = R_specific_air * Trixi.temperature(u, equations)
+
+        C_air = 120.0
+        mu_ref_air = 1.827e-5
+
+        return mu_ref_air * (T_ref + C_air) / (T + C_air) * (T / T_ref)^1.5
+    end
+
+    function mu_control(u, equations, T_ref, R_specific, C, mu_ref)
+        T = R_specific * Trixi.temperature(u, equations)
+
+        return mu_ref * (T_ref + C) / (T + C) * (T / T_ref)^1.5
+    end
+
+    # Dry air (values from Wikipedia: https://de.wikipedia.org/wiki/Sutherland-Modell)
+    T_ref = 291.15
+    C = 120.0 # Sutherland's constant
+    R_specific = 287.052874
+    mu_ref = 1.827e-5
+    prandtl_number() = 0.72
+    gamma = 1.4
+
+    equations = CompressibleEulerEquations2D(gamma)
+    equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, mu = mu,
+                                                              Prandtl = prandtl_number())
+
+    # Flow at rest
+    u = prim2cons(SVector(1.0, 0.0, 0.0, 1.0), equations_parabolic)
+
+    # Comparison value from https://www.engineeringtoolbox.com/air-absolute-kinematic-viscosity-d_601.html at 18°C
+    @test isapprox(mu_control(u, equations_parabolic, T_ref, R_specific, C, mu_ref),
+                   1.803e-5, atol = 5e-8)
+end
 end
 
 end #module