From c910b288d09e04268be700eafff3bca5d7190ad0 Mon Sep 17 00:00:00 2001
From: Daniel_Doehring <daniel.doehring@rwth-aachen.de>
Date: Mon, 8 Dec 2025 10:01:24 +0100
Subject: [PATCH 01/11] Diffusive CFL `P4est`

---
 ...xir_advection_diffusion_nonperiodic_amr.jl |   3 +-
 .../elixir_navierstokes_vortex_street.jl      |   8 +-
 .../elixir_advection_diffusion_nonperiodic.jl |  89 +++++++++++++
 ...elixir_navierstokes_taylor_green_vortex.jl |   8 +-
 src/callbacks_step/stepsize_dg1d.jl           |  40 +++---
 src/callbacks_step/stepsize_dg2d.jl           |  89 ++++++++++++-
 src/callbacks_step/stepsize_dg3d.jl           | 123 ++++++++++++++----
 test/test_parabolic_2d.jl                     |  44 +++++++
 test/test_parabolic_3d.jl                     |  40 ++++++
 9 files changed, 386 insertions(+), 58 deletions(-)
 create mode 100644 examples/p4est_3d_dgsem/elixir_advection_diffusion_nonperiodic.jl

diff --git a/examples/p4est_2d_dgsem/elixir_advection_diffusion_nonperiodic_amr.jl b/examples/p4est_2d_dgsem/elixir_advection_diffusion_nonperiodic_amr.jl
index 35effcf9729..5e68250899e 100644
--- a/examples/p4est_2d_dgsem/elixir_advection_diffusion_nonperiodic_amr.jl
+++ b/examples/p4est_2d_dgsem/elixir_advection_diffusion_nonperiodic_amr.jl
@@ -94,6 +94,7 @@ callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, amr
 # run the simulation
 
 # OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks
+ode_alg = RDPK3SpFSAL49()
 time_int_tol = 1.0e-11
-sol = solve(ode, dt = 1e-7, RDPK3SpFSAL49(); abstol = time_int_tol, reltol = time_int_tol,
+sol = solve(ode, dt = 1e-7, ode_alg; abstol = time_int_tol, reltol = time_int_tol,
             ode_default_options()..., callback = callbacks)
diff --git a/examples/p4est_2d_dgsem/elixir_navierstokes_vortex_street.jl b/examples/p4est_2d_dgsem/elixir_navierstokes_vortex_street.jl
index 3a868f1dac8..4bf9d0c3ec4 100644
--- a/examples/p4est_2d_dgsem/elixir_navierstokes_vortex_street.jl
+++ b/examples/p4est_2d_dgsem/elixir_navierstokes_vortex_street.jl
@@ -157,9 +157,11 @@ callbacks = CallbackSet(summary_callback,
 ###############################################################################
 # run the simulation
 
+# Moderate number of threads (e.g. 4) advisable to speed things up
+ode_alg = RDPK3SpFSAL49(thread = Trixi.True())
 time_int_tol = 1e-7
-sol = solve(ode,
-            # Moderate number of threads (e.g. 4) advisable to speed things up
-            RDPK3SpFSAL49(thread = Trixi.True());
+sol = solve(ode, ode_alg;
+            # not necessary, added for overwriting in tests
+            adaptive = true, dt = 1.4e-3,
             abstol = time_int_tol, reltol = time_int_tol,
             ode_default_options()..., callback = callbacks)
diff --git a/examples/p4est_3d_dgsem/elixir_advection_diffusion_nonperiodic.jl b/examples/p4est_3d_dgsem/elixir_advection_diffusion_nonperiodic.jl
new file mode 100644
index 00000000000..3157ad4e345
--- /dev/null
+++ b/examples/p4est_3d_dgsem/elixir_advection_diffusion_nonperiodic.jl
@@ -0,0 +1,89 @@
+using OrdinaryDiffEqLowStorageRK
+using Trixi
+
+###############################################################################
+# semidiscretization of the linear advection-diffusion equation
+
+diffusivity() = 1 / 17
+advection_velocity = (1.0, 0.0, 0.0)
+equations = LinearScalarAdvectionEquation3D(advection_velocity)
+equations_parabolic = LaplaceDiffusion3D(diffusivity(), equations)
+
+polydeg = 3
+solver = DGSEM(polydeg = polydeg, surface_flux = flux_lax_friedrichs)
+
+coordinates_min = (-1.0, -0.5, -0.5)
+coordinates_max = (0.0, 0.5, 0.5)
+
+# Affine type mapping to take the [-1,1]^3 domain
+# and warp it as described in https://arxiv.org/abs/2012.12040
+function mapping(xi, eta, zeta)
+    y_ = eta + 1 / 6 * (cos(1.5 * pi * xi) * cos(0.5 * pi * eta) * cos(0.5 * pi * zeta))
+    x_ = xi + 1 / 6 * (cos(0.5 * pi * xi) * cos(2 * pi * y_) * cos(0.5 * pi * zeta))
+    z_ = zeta + 1 / 6 * (cos(0.5 * pi * x_) * cos(pi * y_) * cos(0.5 * pi * zeta))
+
+    # Map from [-1, 1]^3 to [-1, -0.5, -0.5] x [0, 0.5, 0.5]
+    x = 0.5 * (x_ - 1)
+    y = 0.5 * y_
+    z = 0.5 * z_
+
+    return SVector(x, y, z)
+end
+
+trees_per_dimension = (5, 3, 3)
+mesh = P4estMesh(trees_per_dimension, polydeg = polydeg,
+                 mapping = mapping, periodicity = false)
+
+# Example setup taken from
+# - Truman Ellis, Jesse Chan, and Leszek Demkowicz (2016).
+#   Robust DPG methods for transient convection-diffusion.
+#   In: Building bridges: connections and challenges in modern approaches
+#   to numerical partial differential equations.
+#   [DOI](https://doi.org/10.1007/978-3-319-41640-3_6).
+function initial_condition_eriksson_johnson(x, t, equations)
+    l = 4
+    epsilon = diffusivity() # NOTE: this requires epsilon <= 1/16 due to sqrt
+    lambda_1 = (-1 + sqrt(1 - 4 * epsilon * l)) / (-2 * epsilon)
+    lambda_2 = (-1 - sqrt(1 - 4 * epsilon * l)) / (-2 * epsilon)
+    r1 = (1 + sqrt(1 + 4 * pi^2 * epsilon^2)) / (2 * epsilon)
+    s1 = (1 - sqrt(1 + 4 * pi^2 * epsilon^2)) / (2 * epsilon)
+    u = exp(-l * t) * (exp(lambda_1 * x[1]) - exp(lambda_2 * x[1])) +
+        cos(pi * x[2]) * (exp(s1 * x[1]) - exp(r1 * x[1])) / (exp(-s1) - exp(-r1))
+    return SVector(u)
+end
+initial_condition = initial_condition_eriksson_johnson
+
+boundary_conditions = boundary_condition_default(mesh,
+                                                 BoundaryConditionDirichlet(initial_condition))
+
+semi = SemidiscretizationHyperbolicParabolic(mesh,
+                                             (equations, equations_parabolic),
+                                             initial_condition, solver;
+                                             boundary_conditions = (boundary_conditions,
+                                                                    boundary_conditions))
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 0.5)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 100
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+stepsize_callback = StepsizeCallback(cfl = 1.6,
+                                     cfl_diffusive = 0.25)
+
+callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback,
+                        stepsize_callback)
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false),
+            dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
+            save_everystep = false, callback = callbacks);
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 42aaf1090b6..00e513c3231 100644
--- a/examples/p4est_3d_dgsem/elixir_navierstokes_taylor_green_vortex.jl
+++ b/examples/p4est_3d_dgsem/elixir_navierstokes_taylor_green_vortex.jl
@@ -69,7 +69,7 @@ save_solution = SaveSolutionCallback(interval = 100,
                                      save_initial_solution = true,
                                      save_final_solution = true,
                                      solution_variables = cons2prim)
-alive_callback = AliveCallback(analysis_interval = analysis_interval)
+alive_callback = AliveCallback(alive_interval = 1)
 
 callbacks = CallbackSet(summary_callback,
                         analysis_callback,
@@ -78,6 +78,10 @@ callbacks = CallbackSet(summary_callback,
 ###############################################################################
 # run the simulation
 
+ode_alg = RDPK3SpFSAL49()
 time_int_tol = 1e-8
-sol = solve(ode, RDPK3SpFSAL49(); abstol = time_int_tol, reltol = time_int_tol,
+sol = solve(ode, ode_alg;
+            # not necessary, added for overwriting in tests
+            adaptive = true, dt = 3.913e-3,
+            abstol = time_int_tol, reltol = time_int_tol,
             ode_default_options()..., callback = callbacks)
diff --git a/src/callbacks_step/stepsize_dg1d.jl b/src/callbacks_step/stepsize_dg1d.jl
index a94e80647d6..61d7557919b 100644
--- a/src/callbacks_step/stepsize_dg1d.jl
+++ b/src/callbacks_step/stepsize_dg1d.jl
@@ -8,8 +8,7 @@
 function max_dt(u, t, mesh::TreeMesh{1},
                 constant_speed::False, equations,
                 dg::DG, cache)
-    # to avoid a division by zero if the speed vanishes everywhere,
-    # e.g. for steady-state linear advection
+    # Avoid division by zero if the speed vanishes everywhere
     max_scaled_speed = nextfloat(zero(t))
 
     @batch reduction=(max, max_scaled_speed) for element in eachelement(dg, cache)
@@ -33,28 +32,31 @@ function max_dt(u, t, mesh::TreeMesh{1},
                 constant_diffusivity::False, equations,
                 equations_parabolic::AbstractEquationsParabolic,
                 dg::DG, cache)
-    # to avoid a division by zero if the diffusivity vanishes everywhere
-    max_scaled_speed = nextfloat(zero(t))
+    # Avoid division by zero if the diffusivity vanishes everywhere
+    max_scaled_diffusivity = nextfloat(zero(t))
 
-    @batch reduction=(max, max_scaled_speed) for element in eachelement(dg, cache)
+    @batch reduction=(max, max_scaled_diffusivity) for element in eachelement(dg, cache)
         max_diffusivity_ = zero(max_scaled_speed)
         for i in eachnode(dg)
             u_node = get_node_vars(u, equations, dg, i, element)
             diffusivity = max_diffusivity(u_node, equations_parabolic)
-            max_diffusivity_ = max(max_diffusivity_, diffusivity)
+            max_diffusivity_ = Base.max(max_diffusivity_, diffusivity)
         end
         inv_jacobian = cache.elements.inverse_jacobian[element] # 2 / Δx
-        max_scaled_speed = max(max_scaled_speed, inv_jacobian^2 * max_diffusivity_)
+        # Use `Base.max` to prevent silent failures, as `max` from `@fastmath` doesn't propagate
+        # `NaN`s properly. See https://github.com/trixi-framework/Trixi.jl/pull/2445#discussion_r2336812323
+        max_scaled_diffusivity = Base.max(max_scaled_diffusivity,
+                                          inv_jacobian^2 * max_diffusivity_)
     end
 
     # Factor 4 cancels with 2^2 from `inv_jacobian^2`, resulting in Δx^2
-    return 4 / (nnodes(dg) * max_scaled_speed)
+    return 4 / (nnodes(dg) * max_scaled_diffusivity)
 end
 
 function max_dt(u, t, mesh::TreeMesh{1},
                 constant_speed::True, equations,
                 dg::DG, cache)
-    # to avoid a division by zero if the speed vanishes everywhere,
+    # Avoid division by zero if the speed vanishes everywhere,
     # e.g. for steady-state linear advection
     max_scaled_speed = nextfloat(zero(t))
 
@@ -71,32 +73,32 @@ function max_dt(u, t, mesh::TreeMesh{1},
     return 2 / (nnodes(dg) * max_scaled_speed)
 end
 
-function max_dt(u, t, mesh::TreeMesh,
+function max_dt(u, t, mesh::TreeMesh, # for all dimensions
                 constant_diffusivity::True, equations,
                 equations_parabolic::AbstractEquationsParabolic,
                 dg::DG, cache)
-    # to avoid a division by zero if the diffusivity vanishes everywhere
-    max_scaled_speed = nextfloat(zero(t))
+    # Avoid division by zero if the diffusivity vanishes everywhere
+    max_scaled_diffusivity = nextfloat(zero(t))
 
     diffusivity = max_diffusivity(equations_parabolic)
 
-    @batch reduction=(max, max_scaled_speed) for element in eachelement(dg, cache)
+    @batch reduction=(max, max_scaled_diffusivity) for element in eachelement(dg, cache)
         inv_jacobian = cache.elements.inverse_jacobian[element] # 2 / Δx
         # Note: For the currently supported parabolic equations
-        # Diffusion & Navier-Stokes, we only have one diffusivity,
+        # Diffusion & Navier-Stokes, we only have one isotropic diffusivity,
         # so this is valid for 1D, 2D and 3D.
-        max_scaled_speed = max(max_scaled_speed, inv_jacobian^2 * diffusivity)
+        max_scaled_diffusivity = Base.max(max_scaled_diffusivity,
+                                          inv_jacobian^2 * diffusivity)
     end
 
     # Factor 4 cancels with 2^2 from `inv_jacobian^2`, resulting in Δx^2
-    return 4 / (nnodes(dg) * max_scaled_speed)
+    return 4 / (nnodes(dg) * max_scaled_diffusivity)
 end
 
 function max_dt(u, t, mesh::StructuredMesh{1},
                 constant_speed::False, equations,
                 dg::DG, cache)
-    # to avoid a division by zero if the speed vanishes everywhere,
-    # e.g. for steady-state linear advection
+    # Avoid division by zero if the speed vanishes everywhere
     max_scaled_speed = nextfloat(zero(t))
 
     @batch reduction=(max, max_scaled_speed) for element in eachelement(dg, cache)
@@ -123,7 +125,7 @@ end
 function max_dt(u, t, mesh::StructuredMesh{1},
                 constant_speed::True, equations,
                 dg::DG, cache)
-    # to avoid a division by zero if the speed vanishes everywhere,
+    # Avoid division by zero if the speed vanishes everywhere,
     # e.g. for steady-state linear advection
     max_scaled_speed = nextfloat(zero(t))
 
diff --git a/src/callbacks_step/stepsize_dg2d.jl b/src/callbacks_step/stepsize_dg2d.jl
index 28268be65bd..4c901196e06 100644
--- a/src/callbacks_step/stepsize_dg2d.jl
+++ b/src/callbacks_step/stepsize_dg2d.jl
@@ -7,7 +7,7 @@
 
 function max_dt(u, t, mesh::TreeMesh{2},
                 constant_speed::False, equations, dg::DG, cache)
-    # to avoid a division by zero if the speed vanishes everywhere,
+    # Avoid division by zero if the speed vanishes everywhere,
     # e.g. for steady-state linear advection
     max_scaled_speed = nextfloat(zero(t))
 
@@ -33,7 +33,7 @@ function max_dt(u, t, mesh::TreeMesh{2},
                 constant_diffusivity::False, equations,
                 equations_parabolic::AbstractEquationsParabolic,
                 dg::DG, cache)
-    # to avoid a division by zero if the diffusivity vanishes everywhere
+    # Avoid division by zero if the diffusivity vanishes everywhere
     max_scaled_speed = nextfloat(zero(t))
 
     @batch reduction=(max, max_scaled_speed) for element in eachelement(dg, cache)
@@ -55,7 +55,7 @@ end
 
 function max_dt(u, t, mesh::TreeMesh{2},
                 constant_speed::True, equations, dg::DG, cache)
-    # to avoid a division by zero if the speed vanishes everywhere,
+    # Avoid division by zero if the speed vanishes everywhere,
     # e.g. for steady-state linear advection
     max_scaled_speed = nextfloat(zero(t))
 
@@ -114,7 +114,7 @@ function max_dt(u, t,
                 mesh::Union{StructuredMesh{2}, UnstructuredMesh2D, P4estMesh{2},
                             T8codeMesh{2}, StructuredMeshView{2}},
                 constant_speed::False, equations, dg::DG, cache)
-    # to avoid a division by zero if the speed vanishes everywhere,
+    # Avoid division by zero if the speed vanishes everywhere,
     # e.g. for steady-state linear advection
     max_scaled_speed = nextfloat(zero(t))
 
@@ -148,13 +148,53 @@ function max_dt(u, t,
     return 2 / (nnodes(dg) * max_scaled_speed)
 end
 
+function max_dt(u, t,
+                mesh::P4estMesh{2}, # Parabolic terms currently only for `TreeMesh` and `P4estMesh`
+                constant_diffusivity::False, equations,
+                equations_parabolic::AbstractEquationsParabolic,
+                dg::DG, cache)
+    # Avoid division by zero if the diffusivity vanishes everywhere
+    max_scaled_diffusivity = nextfloat(zero(t))
+
+    @unpack contravariant_vectors, inverse_jacobian = cache.elements
+
+    @batch reduction=(max, max_scaled_diffusivity) for element in eachelement(dg, cache)
+        max_diffusivity1 = max_diffusivity2 = zero(max_scaled_diffusivity)
+        for j in eachnode(dg), i in eachnode(dg)
+            u_node = get_node_vars(u, equations, dg, i, j, element)
+            diffusivity = max_diffusivity(u_node, equations_parabolic)
+
+            # Local diffusivity transformed to the reference element
+            Ja11, Ja12 = get_contravariant_vector(1, contravariant_vectors,
+                                                  i, j, element)
+            diffusivity1_transformed = diffusivity * (Ja11^2 + Ja12^2)
+            Ja21, Ja22 = get_contravariant_vector(2, contravariant_vectors,
+                                                  i, j, element)
+            diffusivity2_transformed = diffusivity * (Ja21^2 + Ja22^2)
+
+            inv_jacobian = abs(inverse_jacobian[i, j, element])
+
+            max_diffusivity1 = Base.max(max_diffusivity1,
+                                        diffusivity1_transformed * inv_jacobian^2)
+            max_diffusivity2 = Base.max(max_diffusivity2,
+                                        diffusivity2_transformed * inv_jacobian^2)
+        end
+        # Use `Base.max` to prevent silent failures, as `max` from `@fastmath` doesn't propagate
+        # `NaN`s properly. See https://github.com/trixi-framework/Trixi.jl/pull/2445#discussion_r2336812323
+        max_scaled_diffusivity = Base.max(max_scaled_diffusivity,
+                                          max_diffusivity1 + max_diffusivity2)
+    end
+
+    return 4 / (nnodes(dg) * max_scaled_diffusivity)
+end
+
 function max_dt(u, t,
                 mesh::Union{StructuredMesh{2}, UnstructuredMesh2D, P4estMesh{2},
                             P4estMeshView{2}, T8codeMesh{2}, StructuredMeshView{2}},
                 constant_speed::True, equations, dg::DG, cache)
     @unpack contravariant_vectors, inverse_jacobian = cache.elements
 
-    # to avoid a division by zero if the speed vanishes everywhere,
+    # Avoid division by zero if the speed vanishes everywhere,
     # e.g. for steady-state linear advection
     max_scaled_speed = nextfloat(zero(t))
 
@@ -182,6 +222,45 @@ function max_dt(u, t,
     return 2 / (nnodes(dg) * max_scaled_speed)
 end
 
+function max_dt(u, t,
+                mesh::P4estMesh{2}, # Parabolic terms currently only for `TreeMesh` and `P4estMesh`
+                constant_diffusivity::True, equations,
+                equations_parabolic::AbstractEquationsParabolic,
+                dg::DG, cache)
+    # Avoid division by zero if the diffusivity vanishes everywhere
+    max_scaled_diffusivity = nextfloat(zero(t))
+
+    diffusivity = max_diffusivity(equations_parabolic)
+
+    @unpack contravariant_vectors, inverse_jacobian = cache.elements
+
+    @batch reduction=(max, max_scaled_diffusivity) for element in eachelement(dg, cache)
+        max_diffusivity1 = max_diffusivity2 = zero(max_scaled_diffusivity)
+        for j in eachnode(dg), i in eachnode(dg)
+            # Local diffusivity transformed to the reference element
+            Ja11, Ja12 = get_contravariant_vector(1, contravariant_vectors,
+                                                  i, j, element)
+            diffusivity1_transformed = diffusivity * (Ja11^2 + Ja12^2)
+            Ja21, Ja22 = get_contravariant_vector(2, contravariant_vectors,
+                                                  i, j, element)
+            diffusivity2_transformed = diffusivity * (Ja21^2 + Ja22^2)
+
+            inv_jacobian = abs(inverse_jacobian[i, j, element])
+
+            max_diffusivity1 = Base.max(max_diffusivity1,
+                                        diffusivity1_transformed * inv_jacobian^2)
+            max_diffusivity2 = Base.max(max_diffusivity2,
+                                        diffusivity2_transformed * inv_jacobian^2)
+        end
+        # Use `Base.max` to prevent silent failures, as `max` from `@fastmath` doesn't propagate
+        # `NaN`s properly. See https://github.com/trixi-framework/Trixi.jl/pull/2445#discussion_r2336812323
+        max_scaled_diffusivity = Base.max(max_scaled_diffusivity,
+                                          max_diffusivity1 + max_diffusivity2)
+    end
+
+    return 4 / (nnodes(dg) * max_scaled_diffusivity)
+end
+
 function max_dt(u, t, mesh::ParallelP4estMesh{2},
                 constant_speed::False, equations, dg::DG, cache)
     # call the method accepting a general `mesh::P4estMesh{2}`
diff --git a/src/callbacks_step/stepsize_dg3d.jl b/src/callbacks_step/stepsize_dg3d.jl
index dba95999d3a..61fe15e13b3 100644
--- a/src/callbacks_step/stepsize_dg3d.jl
+++ b/src/callbacks_step/stepsize_dg3d.jl
@@ -7,7 +7,7 @@
 
 function max_dt(u, t, mesh::TreeMesh{3},
                 constant_speed::False, equations, dg::DG, cache)
-    # to avoid a division by zero if the speed vanishes everywhere,
+    # Avoid division by zero if the speed vanishes everywhere,
     # e.g. for steady-state linear advection
     max_scaled_speed = nextfloat(zero(t))
 
@@ -31,33 +31,9 @@ function max_dt(u, t, mesh::TreeMesh{3},
     return 2 / (nnodes(dg) * max_scaled_speed)
 end
 
-function max_dt(u, t, mesh::TreeMesh{3},
-                constant_diffusivity::False, equations,
-                equations_parabolic::AbstractEquationsParabolic,
-                dg::DG, cache)
-    # to avoid a division by zero if the diffusivity vanishes everywhere
-    max_scaled_speed = nextfloat(zero(t))
-
-    @batch reduction=(max, max_scaled_speed) for element in eachelement(dg, cache)
-        max_diffusivity_ = zero(max_scaled_speed)
-        for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg)
-            u_node = get_node_vars(u, equations, dg, i, j, k, element)
-            # Note: For the currently supported parabolic equations
-            # Diffusion & Navier-Stokes, we only have one diffusivity.
-            diffusivity = max_diffusivity(u_node, equations_parabolic)
-            max_diffusivity_ = max(max_diffusivity_, diffusivity)
-        end
-        inv_jacobian = cache.elements.inverse_jacobian[element] # 2 / Δx
-        max_scaled_speed = max(max_scaled_speed, inv_jacobian^2 * max_diffusivity_)
-    end
-
-    # Factor 4 cancels with 2^2 from `inv_jacobian^2`, resulting in Δx^2
-    return 4 / (nnodes(dg) * max_scaled_speed)
-end
-
 function max_dt(u, t, mesh::TreeMesh{3},
                 constant_speed::True, equations, dg::DG, cache)
-    # to avoid a division by zero if the speed vanishes everywhere,
+    # Avoid division by zero if the speed vanishes everywhere,
     # e.g. for steady-state linear advection
     max_scaled_speed = nextfloat(zero(t))
 
@@ -77,7 +53,7 @@ end
 
 function max_dt(u, t, mesh::Union{StructuredMesh{3}, P4estMesh{3}, T8codeMesh{3}},
                 constant_speed::False, equations, dg::DG, cache)
-    # to avoid a division by zero if the speed vanishes everywhere,
+    # Avoid division by zero if the speed vanishes everywhere,
     # e.g. for steady-state linear advection
     max_scaled_speed = nextfloat(zero(t))
 
@@ -115,9 +91,55 @@ function max_dt(u, t, mesh::Union{StructuredMesh{3}, P4estMesh{3}, T8codeMesh{3}
     return 2 / (nnodes(dg) * max_scaled_speed)
 end
 
+function max_dt(u, t,
+                mesh::P4estMesh{3}, # Parabolic terms currently only for `TreeMesh` and `P4estMesh`
+                constant_diffusivity::False, equations,
+                equations_parabolic::AbstractEquationsParabolic,
+                dg::DG, cache)
+    # Avoid division by zero if the diffusivity vanishes everywhere
+    max_scaled_diffusivity = nextfloat(zero(t))
+
+    @unpack contravariant_vectors, inverse_jacobian = cache.elements
+
+    @batch reduction=(max, max_scaled_diffusivity) for element in eachelement(dg, cache)
+        max_diffusivity1 = max_diffusivity2 = max_diffusivity3 = zero(max_scaled_diffusivity)
+        for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg)
+            u_node = get_node_vars(u, equations, dg, i, j, k, element)
+            diffusivity = max_diffusivity(u_node, equations_parabolic)
+
+            # Local diffusivity transformed to the reference element
+            Ja11, Ja12, Ja13 = get_contravariant_vector(1, contravariant_vectors,
+                                                        i, j, k, element)
+            diffusivity1_transformed = diffusivity * (Ja11^2 + Ja12^2 + Ja13^2)
+            Ja21, Ja22, Ja23 = get_contravariant_vector(2, contravariant_vectors,
+                                                        i, j, k, element)
+            diffusivity2_transformed = diffusivity * (Ja21^2 + Ja22^2 + Ja23^2)
+            Ja31, Ja32, Ja33 = get_contravariant_vector(3, contravariant_vectors,
+                                                        i, j, k, element)
+            diffusivity3_transformed = diffusivity * (Ja31^2 + Ja32^2 + Ja33^2)
+
+            inv_jacobian = abs(inverse_jacobian[i, j, k, element])
+
+            max_diffusivity1 = Base.max(max_diffusivity1,
+                                        diffusivity1_transformed * inv_jacobian^2)
+            max_diffusivity2 = Base.max(max_diffusivity2,
+                                        diffusivity2_transformed * inv_jacobian^2)
+            max_diffusivity3 = Base.max(max_diffusivity3,
+                                        diffusivity3_transformed * inv_jacobian^2)
+        end
+        # Use `Base.max` to prevent silent failures, as `max` from `@fastmath` doesn't propagate
+        # `NaN`s properly. See https://github.com/trixi-framework/Trixi.jl/pull/2445#discussion_r2336812323
+        max_scaled_diffusivity = Base.max(max_scaled_diffusivity,
+                                          max_diffusivity1 + max_diffusivity2 +
+                                          max_diffusivity3)
+    end
+
+    return 4 / (nnodes(dg) * max_scaled_diffusivity)
+end
+
 function max_dt(u, t, mesh::Union{StructuredMesh{3}, P4estMesh{3}, T8codeMesh{3}},
                 constant_speed::True, equations, dg::DG, cache)
-    # to avoid a division by zero if the speed vanishes everywhere,
+    # Avoid division by zero if the speed vanishes everywhere,
     # e.g. for steady-state linear advection
     max_scaled_speed = nextfloat(zero(t))
 
@@ -154,6 +176,51 @@ function max_dt(u, t, mesh::Union{StructuredMesh{3}, P4estMesh{3}, T8codeMesh{3}
     return 2 / (nnodes(dg) * max_scaled_speed)
 end
 
+function max_dt(u, t,
+                mesh::P4estMesh{3}, # Parabolic terms currently only for `TreeMesh` and `P4estMesh`
+                constant_diffusivity::True, equations,
+                equations_parabolic::AbstractEquationsParabolic,
+                dg::DG, cache)
+    # Avoid division by zero if the diffusivity vanishes everywhere
+    max_scaled_diffusivity = nextfloat(zero(t))
+
+    diffusivity = max_diffusivity(equations_parabolic)
+
+    @unpack contravariant_vectors, inverse_jacobian = cache.elements
+
+    @batch reduction=(max, max_scaled_diffusivity) for element in eachelement(dg, cache)
+        max_diffusivity1 = max_diffusivity2 = max_diffusivity3 = zero(max_scaled_diffusivity)
+        for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg)
+            # Local diffusivity transformed to the reference element
+            Ja11, Ja12, Ja13 = get_contravariant_vector(1, contravariant_vectors,
+                                                        i, j, k, element)
+            diffusivity1_transformed = diffusivity * (Ja11^2 + Ja12^2 + Ja13^2)
+            Ja21, Ja22, Ja23 = get_contravariant_vector(2, contravariant_vectors,
+                                                        i, j, k, element)
+            diffusivity2_transformed = diffusivity * (Ja21^2 + Ja22^2 + Ja23^2)
+            Ja31, Ja32, Ja33 = get_contravariant_vector(3, contravariant_vectors,
+                                                        i, j, k, element)
+            diffusivity3_transformed = diffusivity * (Ja31^2 + Ja32^2 + Ja33^2)
+
+            inv_jacobian = abs(inverse_jacobian[i, j, k, element])
+
+            max_diffusivity1 = Base.max(max_diffusivity1,
+                                        diffusivity1_transformed * inv_jacobian^2)
+            max_diffusivity2 = Base.max(max_diffusivity2,
+                                        diffusivity2_transformed * inv_jacobian^2)
+            max_diffusivity3 = Base.max(max_diffusivity3,
+                                        diffusivity3_transformed * inv_jacobian^2)
+        end
+        # Use `Base.max` to prevent silent failures, as `max` from `@fastmath` doesn't propagate
+        # `NaN`s properly. See https://github.com/trixi-framework/Trixi.jl/pull/2445#discussion_r2336812323
+        max_scaled_diffusivity = Base.max(max_scaled_diffusivity,
+                                          max_diffusivity1 + max_diffusivity2 +
+                                          max_diffusivity3)
+    end
+
+    return 4 / (nnodes(dg) * max_scaled_diffusivity)
+end
+
 function max_dt(u, t, mesh::ParallelP4estMesh{3},
                 constant_speed::False, equations, dg::DG, cache)
     # call the method accepting a general `mesh::P4estMesh{3}`
diff --git a/test/test_parabolic_2d.jl b/test/test_parabolic_2d.jl
index f69dcd7df5d..40f395d5956 100644
--- a/test/test_parabolic_2d.jl
+++ b/test/test_parabolic_2d.jl
@@ -565,6 +565,22 @@ end
     @test_allocations(Trixi.rhs!, semi, sol, 1000)
 end
 
+@trixi_testset "P4estMesh2D: elixir_advection_diffusion_nonperiodic_amr.jl (Diffusive CFL)" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR, "p4est_2d_dgsem",
+                                 "elixir_advection_diffusion_nonperiodic_amr.jl"),
+                        callbacks=CallbackSet(summary_callback, analysis_callback,
+                                              alive_callback,
+                                              StepsizeCallback(cfl = 1.6,
+                                                               cfl_diffusive = 0.2)),
+                        ode_alg=CarpenterKennedy2N54(williamson_condition = false),
+                        dt=1.0, # will be overwritten
+                        l2=[0.00010850375815619432],
+                        linf=[0.0024081141187764932])
+    # Ensure that we do not have excessive memory allocations
+    # (e.g., from type instabilities)
+    @test_allocations(Trixi.rhs!, semi, sol, 1000)
+end
+
 @trixi_testset "P4estMesh2D: elixir_advection_diffusion_nonperiodic_curved.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "p4est_2d_dgsem",
                                  "elixir_advection_diffusion_nonperiodic_curved.jl"),
@@ -806,6 +822,34 @@ end
     @test_allocations(Trixi.rhs!, semi, sol, 1000)
 end
 
+@trixi_testset "elixir_navierstokes_vortex_street.jl" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR, "p4est_2d_dgsem",
+                                 "elixir_navierstokes_vortex_street.jl"),
+                        tspan=(0.0, 0.5),
+                        Re=20, # Render flow diffusion-dominated
+                        callbacks=CallbackSet(summary_callback, analysis_callback,
+                                              alive_callback,
+                                              StepsizeCallback(cfl = 2.3,
+                                                               cfl_diffusive = 1.0)),
+                        adaptive=false, # respect CFL
+                        ode_alg=CKLLSRK95_4S(),
+                        l2=[
+                            0.011916725799140692,
+                            0.027926098816747836,
+                            0.01902700347912797,
+                            0.11793406377747188
+                        ],
+                        linf=[
+                            0.3546113252441576,
+                            1.0152021857472098,
+                            0.5811488174143082,
+                            3.207373092525428
+                        ])
+    # Ensure that we do not have excessive memory allocations
+    # (e.g., from type instabilities)
+    @test_allocations(Trixi.rhs!, semi, sol, 1000)
+end
+
 @trixi_testset "elixir_navierstokes_vortex_street.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "p4est_2d_dgsem",
                                  "elixir_navierstokes_vortex_street.jl"),
diff --git a/test/test_parabolic_3d.jl b/test/test_parabolic_3d.jl
index a346671127b..60b2f9a9a9b 100644
--- a/test/test_parabolic_3d.jl
+++ b/test/test_parabolic_3d.jl
@@ -313,6 +313,15 @@ end
     @test_allocations(Trixi.Trixi.rhs_parabolic!, semi, sol, 100)
 end
 
+@trixi_testset "P4estMesh3D: elixir_advection_diffusion_nonperiodic.jl" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR, "p4est_3d_dgsem",
+                                 "elixir_advection_diffusion_nonperiodic.jl"),
+                        l2=[0.006421164728264022], linf=[0.41638021060047015])
+    # Ensure that we do not have excessive memory allocations
+    # (e.g., from type instabilities)
+    @test_allocations(Trixi.rhs!, semi, sol, 1000)
+end
+
 @trixi_testset "P4estMesh3D: elixir_navierstokes_convergence.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "p4est_3d_dgsem",
                                  "elixir_navierstokes_convergence.jl"),
@@ -360,6 +369,37 @@ end
     @test_allocations(Trixi.rhs!, semi, sol, 1000)
 end
 
+@trixi_testset "P4estMesh3D: elixir_navierstokes_taylor_green_vortex.jl" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR, "p4est_3d_dgsem",
+                                 "elixir_navierstokes_taylor_green_vortex.jl"),
+                        initial_refinement_level=2, tspan=(0.0, 0.25),
+                        surface_flux=FluxHLL(min_max_speed_naive),
+                        mu=0.5, # render flow diffusion-dominated
+                        callbacks=CallbackSet(summary_callback, analysis_callback,
+                                              alive_callback,
+                                              StepsizeCallback(cfl = 2.3,
+                                                               cfl_diffusive = 0.4)),
+                        adaptive=false, # respect CFL
+                        ode_alg=CKLLSRK95_4S(),
+                        l2=[
+                            0.00013361138825044455,
+                            0.1110290325891525,
+                            0.11102903258915098,
+                            0.009310736674102591,
+                            0.18713598964308234
+                        ],
+                        linf=[
+                            0.0005501488385279973,
+                            0.3151241675791191,
+                            0.3151241709442578,
+                            0.026332750964829052,
+                            0.5712639505309198
+                        ])
+    # Ensure that we do not have excessive memory allocations
+    # (e.g., from type instabilities)
+    @test_allocations(Trixi.rhs!, semi, sol, 1000)
+end
+
 @trixi_testset "TreeMesh3D: elixir_advection_diffusion_amr.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "tree_3d_dgsem",
                                  "elixir_advection_diffusion_amr.jl"),

From 118afe81e1cbf4b2f6cc88f5dadf0d34adf2e5a4 Mon Sep 17 00:00:00 2001
From: Daniel_Doehring <daniel.doehring@rwth-aachen.de>
Date: Mon, 8 Dec 2025 10:03:18 +0100
Subject: [PATCH 02/11] rev

---
 .../p4est_3d_dgsem/elixir_navierstokes_taylor_green_vortex.jl   | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

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 00e513c3231..8a8dc1f8582 100644
--- a/examples/p4est_3d_dgsem/elixir_navierstokes_taylor_green_vortex.jl
+++ b/examples/p4est_3d_dgsem/elixir_navierstokes_taylor_green_vortex.jl
@@ -69,7 +69,7 @@ save_solution = SaveSolutionCallback(interval = 100,
                                      save_initial_solution = true,
                                      save_final_solution = true,
                                      solution_variables = cons2prim)
-alive_callback = AliveCallback(alive_interval = 1)
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
 
 callbacks = CallbackSet(summary_callback,
                         analysis_callback,

From 453cffcb378ade71b808f10d9951a1ea042e6855 Mon Sep 17 00:00:00 2001
From: Daniel_Doehring <daniel.doehring@rwth-aachen.de>
Date: Mon, 8 Dec 2025 10:10:26 +0100
Subject: [PATCH 03/11] rev

---
 src/callbacks_step/stepsize_dg3d.jl | 24 ++++++++++++++++++++++++
 1 file changed, 24 insertions(+)

diff --git a/src/callbacks_step/stepsize_dg3d.jl b/src/callbacks_step/stepsize_dg3d.jl
index 61fe15e13b3..f242004073e 100644
--- a/src/callbacks_step/stepsize_dg3d.jl
+++ b/src/callbacks_step/stepsize_dg3d.jl
@@ -51,6 +51,30 @@ function max_dt(u, t, mesh::TreeMesh{3},
     return 2 / (nnodes(dg) * max_scaled_speed)
 end
 
+function max_dt(u, t, mesh::TreeMesh{3},
+                constant_diffusivity::False, equations,
+                equations_parabolic::AbstractEquationsParabolic,
+                dg::DG, cache)
+    # to avoid a division by zero if the diffusivity vanishes everywhere
+    max_scaled_speed = nextfloat(zero(t))
+
+    @batch reduction=(max, max_scaled_speed) for element in eachelement(dg, cache)
+        max_diffusivity_ = zero(max_scaled_speed)
+        for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg)
+            u_node = get_node_vars(u, equations, dg, i, j, k, element)
+            # Note: For the currently supported parabolic equations
+            # Diffusion & Navier-Stokes, we only have one diffusivity.
+            diffusivity = max_diffusivity(u_node, equations_parabolic)
+            max_diffusivity_ = max(max_diffusivity_, diffusivity)
+        end
+        inv_jacobian = cache.elements.inverse_jacobian[element] # 2 / Δx
+        max_scaled_speed = max(max_scaled_speed, inv_jacobian^2 * max_diffusivity_)
+    end
+
+    # Factor 4 cancels with 2^2 from `inv_jacobian^2`, resulting in Δx^2
+    return 4 / (nnodes(dg) * max_scaled_speed)
+end
+
 function max_dt(u, t, mesh::Union{StructuredMesh{3}, P4estMesh{3}, T8codeMesh{3}},
                 constant_speed::False, equations, dg::DG, cache)
     # Avoid division by zero if the speed vanishes everywhere,

From 7135f2e66530d0ee708070aa864e23bfe9d44052 Mon Sep 17 00:00:00 2001
From: Daniel_Doehring <daniel.doehring@rwth-aachen.de>
Date: Mon, 8 Dec 2025 10:11:13 +0100
Subject: [PATCH 04/11] up

---
 src/callbacks_step/stepsize_dg3d.jl | 13 +++++++------
 1 file changed, 7 insertions(+), 6 deletions(-)

diff --git a/src/callbacks_step/stepsize_dg3d.jl b/src/callbacks_step/stepsize_dg3d.jl
index f242004073e..8f2bc43efa4 100644
--- a/src/callbacks_step/stepsize_dg3d.jl
+++ b/src/callbacks_step/stepsize_dg3d.jl
@@ -55,11 +55,11 @@ function max_dt(u, t, mesh::TreeMesh{3},
                 constant_diffusivity::False, equations,
                 equations_parabolic::AbstractEquationsParabolic,
                 dg::DG, cache)
-    # to avoid a division by zero if the diffusivity vanishes everywhere
-    max_scaled_speed = nextfloat(zero(t))
+    # Avoid division by zero if the diffusivity vanishes everywhere
+    max_scaled_diffusivity = nextfloat(zero(t))
 
-    @batch reduction=(max, max_scaled_speed) for element in eachelement(dg, cache)
-        max_diffusivity_ = zero(max_scaled_speed)
+    @batch reduction=(max, max_scaled_diffusivity) for element in eachelement(dg, cache)
+        max_diffusivity_ = zero(max_scaled_diffusivity)
         for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg)
             u_node = get_node_vars(u, equations, dg, i, j, k, element)
             # Note: For the currently supported parabolic equations
@@ -68,11 +68,12 @@ function max_dt(u, t, mesh::TreeMesh{3},
             max_diffusivity_ = max(max_diffusivity_, diffusivity)
         end
         inv_jacobian = cache.elements.inverse_jacobian[element] # 2 / Δx
-        max_scaled_speed = max(max_scaled_speed, inv_jacobian^2 * max_diffusivity_)
+        max_scaled_diffusivity = max(max_scaled_diffusivity,
+                                     inv_jacobian^2 * max_diffusivity_)
     end
 
     # Factor 4 cancels with 2^2 from `inv_jacobian^2`, resulting in Δx^2
-    return 4 / (nnodes(dg) * max_scaled_speed)
+    return 4 / (nnodes(dg) * max_scaled_diffusivity)
 end
 
 function max_dt(u, t, mesh::Union{StructuredMesh{3}, P4estMesh{3}, T8codeMesh{3}},

From 5b4693180996da009666765851dce314d4eccdc9 Mon Sep 17 00:00:00 2001
From: Daniel_Doehring <daniel.doehring@rwth-aachen.de>
Date: Mon, 8 Dec 2025 10:13:20 +0100
Subject: [PATCH 05/11] move

---
 src/callbacks_step/stepsize_dg3d.jl | 40 ++++++++++++++---------------
 1 file changed, 20 insertions(+), 20 deletions(-)

diff --git a/src/callbacks_step/stepsize_dg3d.jl b/src/callbacks_step/stepsize_dg3d.jl
index 8f2bc43efa4..107ec8c3e14 100644
--- a/src/callbacks_step/stepsize_dg3d.jl
+++ b/src/callbacks_step/stepsize_dg3d.jl
@@ -31,26 +31,6 @@ function max_dt(u, t, mesh::TreeMesh{3},
     return 2 / (nnodes(dg) * max_scaled_speed)
 end
 
-function max_dt(u, t, mesh::TreeMesh{3},
-                constant_speed::True, equations, dg::DG, cache)
-    # Avoid division by zero if the speed vanishes everywhere,
-    # e.g. for steady-state linear advection
-    max_scaled_speed = nextfloat(zero(t))
-
-    max_lambda1, max_lambda2, max_lambda3 = max_abs_speeds(equations)
-
-    @batch reduction=(max, max_scaled_speed) for element in eachelement(dg, cache)
-        inv_jacobian = cache.elements.inverse_jacobian[element]
-        # Use `Base.max` to prevent silent failures, as `max` from `@fastmath` doesn't propagate
-        # `NaN`s properly. See https://github.com/trixi-framework/Trixi.jl/pull/2445#discussion_r2336812323
-        max_scaled_speed = Base.max(max_scaled_speed,
-                                    inv_jacobian *
-                                    (max_lambda1 + max_lambda2 + max_lambda3))
-    end
-
-    return 2 / (nnodes(dg) * max_scaled_speed)
-end
-
 function max_dt(u, t, mesh::TreeMesh{3},
                 constant_diffusivity::False, equations,
                 equations_parabolic::AbstractEquationsParabolic,
@@ -76,6 +56,26 @@ function max_dt(u, t, mesh::TreeMesh{3},
     return 4 / (nnodes(dg) * max_scaled_diffusivity)
 end
 
+function max_dt(u, t, mesh::TreeMesh{3},
+                constant_speed::True, equations, dg::DG, cache)
+    # Avoid division by zero if the speed vanishes everywhere,
+    # e.g. for steady-state linear advection
+    max_scaled_speed = nextfloat(zero(t))
+
+    max_lambda1, max_lambda2, max_lambda3 = max_abs_speeds(equations)
+
+    @batch reduction=(max, max_scaled_speed) for element in eachelement(dg, cache)
+        inv_jacobian = cache.elements.inverse_jacobian[element]
+        # Use `Base.max` to prevent silent failures, as `max` from `@fastmath` doesn't propagate
+        # `NaN`s properly. See https://github.com/trixi-framework/Trixi.jl/pull/2445#discussion_r2336812323
+        max_scaled_speed = Base.max(max_scaled_speed,
+                                    inv_jacobian *
+                                    (max_lambda1 + max_lambda2 + max_lambda3))
+    end
+
+    return 2 / (nnodes(dg) * max_scaled_speed)
+end
+
 function max_dt(u, t, mesh::Union{StructuredMesh{3}, P4estMesh{3}, T8codeMesh{3}},
                 constant_speed::False, equations, dg::DG, cache)
     # Avoid division by zero if the speed vanishes everywhere,

From be664dd9e00bda6f2ebb9124b03ece6b7b4b520c Mon Sep 17 00:00:00 2001
From: Daniel Doehring <doehringd2@gmail.com>
Date: Mon, 8 Dec 2025 10:49:44 +0100
Subject: [PATCH 06/11] Apply suggestions from code review

---
 test/test_parabolic_3d.jl | 16 ++++++++--------
 1 file changed, 8 insertions(+), 8 deletions(-)

diff --git a/test/test_parabolic_3d.jl b/test/test_parabolic_3d.jl
index 60b2f9a9a9b..05130ebe067 100644
--- a/test/test_parabolic_3d.jl
+++ b/test/test_parabolic_3d.jl
@@ -369,7 +369,7 @@ end
     @test_allocations(Trixi.rhs!, semi, sol, 1000)
 end
 
-@trixi_testset "P4estMesh3D: elixir_navierstokes_taylor_green_vortex.jl" begin
+@trixi_testset "P4estMesh3D: elixir_navierstokes_taylor_green_vortex.jl (Diffusive CFL)" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "p4est_3d_dgsem",
                                  "elixir_navierstokes_taylor_green_vortex.jl"),
                         initial_refinement_level=2, tspan=(0.0, 0.25),
@@ -382,18 +382,18 @@ end
                         adaptive=false, # respect CFL
                         ode_alg=CKLLSRK95_4S(),
                         l2=[
-                            0.00013361138825044455,
+                            0.0001336113719363628,
                             0.1110290325891525,
                             0.11102903258915098,
                             0.009310736674102591,
-                            0.18713598964308234
+                            0.18713598666809547
                         ],
                         linf=[
-                            0.0005501488385279973,
-                            0.3151241675791191,
-                            0.3151241709442578,
-                            0.026332750964829052,
-                            0.5712639505309198
+                            0.0005497560773632948,
+                            0.3151224994649052,
+                            0.3151225081914104,
+                            0.02632750316461941,
+                            0.570941570440084
                         ])
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)

From fd344ea1f9526afacd952837ec37bb4f21188bc6 Mon Sep 17 00:00:00 2001
From: Daniel_Doehring <daniel.doehring@rwth-aachen.de>
Date: Mon, 8 Dec 2025 14:00:08 +0100
Subject: [PATCH 07/11] bf

---
 src/callbacks_step/stepsize_dg1d.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/callbacks_step/stepsize_dg1d.jl b/src/callbacks_step/stepsize_dg1d.jl
index 61d7557919b..8a029543575 100644
--- a/src/callbacks_step/stepsize_dg1d.jl
+++ b/src/callbacks_step/stepsize_dg1d.jl
@@ -36,7 +36,7 @@ function max_dt(u, t, mesh::TreeMesh{1},
     max_scaled_diffusivity = nextfloat(zero(t))
 
     @batch reduction=(max, max_scaled_diffusivity) for element in eachelement(dg, cache)
-        max_diffusivity_ = zero(max_scaled_speed)
+        max_diffusivity_ = zero(max_scaled_diffusivity)
         for i in eachnode(dg)
             u_node = get_node_vars(u, equations, dg, i, element)
             diffusivity = max_diffusivity(u_node, equations_parabolic)

From 7d8d1f893e9043fad18df8c5ac213f7a27ae68d7 Mon Sep 17 00:00:00 2001
From: Daniel Doehring <doehringd2@gmail.com>
Date: Mon, 8 Dec 2025 15:06:52 +0100
Subject: [PATCH 08/11] Update test/test_parabolic_3d.jl

---
 test/test_parabolic_3d.jl | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/test/test_parabolic_3d.jl b/test/test_parabolic_3d.jl
index 05130ebe067..c5cc0885e18 100644
--- a/test/test_parabolic_3d.jl
+++ b/test/test_parabolic_3d.jl
@@ -389,11 +389,11 @@ end
                             0.18713598666809547
                         ],
                         linf=[
-                            0.0005497560773632948,
-                            0.3151224994649052,
-                            0.3151225081914104,
-                            0.02632750316461941,
-                            0.570941570440084
+                            0.0005497841741080034,
+                            0.31512261528156615,
+                            0.31512262804742686,
+                            0.02632785741783134,
+                            0.5709646222385913
                         ])
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)

From c4855d15ce16bd104aadb3b4dda0690895d4fe83 Mon Sep 17 00:00:00 2001
From: Daniel_Doehring <daniel.doehring@rwth-aachen.de>
Date: Wed, 10 Dec 2025 09:53:54 +0100
Subject: [PATCH 09/11] comments

---
 src/callbacks_step/stepsize_dg2d.jl | 12 ++++++++++++
 src/callbacks_step/stepsize_dg3d.jl | 12 ++++++++++++
 2 files changed, 24 insertions(+)

diff --git a/src/callbacks_step/stepsize_dg2d.jl b/src/callbacks_step/stepsize_dg2d.jl
index 4c901196e06..f56912948c8 100644
--- a/src/callbacks_step/stepsize_dg2d.jl
+++ b/src/callbacks_step/stepsize_dg2d.jl
@@ -167,6 +167,12 @@ function max_dt(u, t,
             # Local diffusivity transformed to the reference element
             Ja11, Ja12 = get_contravariant_vector(1, contravariant_vectors,
                                                   i, j, element)
+            # The metric terms are squared due to the second derivatives in the parabolic terms,
+            # which lead to application of the chain rule (coordinate transformation) twice.
+            # See for instance also the implementation in FLEXI
+            # https://github.com/flexi-framework/flexi/blob/e980e8635e6605daf906fc176c9897314a9cd9b1/src/equations/navierstokes/calctimestep.f90#L75-L77
+            # or FLUXO:
+            # https://github.com/project-fluxo/fluxo/blob/c7e0cc9b7fd4569dcab67bbb6e5a25c0a84859f1/src/equation/navierstokes/calctimestep.f90#L130-L132
             diffusivity1_transformed = diffusivity * (Ja11^2 + Ja12^2)
             Ja21, Ja22 = get_contravariant_vector(2, contravariant_vectors,
                                                   i, j, element)
@@ -240,6 +246,12 @@ function max_dt(u, t,
             # Local diffusivity transformed to the reference element
             Ja11, Ja12 = get_contravariant_vector(1, contravariant_vectors,
                                                   i, j, element)
+            # The metric terms are squared due to the second derivatives in the parabolic terms,
+            # which lead to application of the chain rule (coordinate transformation) twice.
+            # See for instance also the implementation in FLEXI
+            # https://github.com/flexi-framework/flexi/blob/e980e8635e6605daf906fc176c9897314a9cd9b1/src/equations/navierstokes/calctimestep.f90#L75-L77
+            # or FLUXO:
+            # https://github.com/project-fluxo/fluxo/blob/c7e0cc9b7fd4569dcab67bbb6e5a25c0a84859f1/src/equation/navierstokes/calctimestep.f90#L130-L132
             diffusivity1_transformed = diffusivity * (Ja11^2 + Ja12^2)
             Ja21, Ja22 = get_contravariant_vector(2, contravariant_vectors,
                                                   i, j, element)
diff --git a/src/callbacks_step/stepsize_dg3d.jl b/src/callbacks_step/stepsize_dg3d.jl
index 107ec8c3e14..47560810782 100644
--- a/src/callbacks_step/stepsize_dg3d.jl
+++ b/src/callbacks_step/stepsize_dg3d.jl
@@ -135,6 +135,12 @@ function max_dt(u, t,
             # Local diffusivity transformed to the reference element
             Ja11, Ja12, Ja13 = get_contravariant_vector(1, contravariant_vectors,
                                                         i, j, k, element)
+            # The metric terms are squared due to the second derivatives in the parabolic terms,
+            # which lead to application of the chain rule (coordinate transformation) twice.
+            # See for instance also the implementation in FLEXI
+            # https://github.com/flexi-framework/flexi/blob/e980e8635e6605daf906fc176c9897314a9cd9b1/src/equations/navierstokes/calctimestep.f90#L75-L77
+            # or FLUXO:
+            # https://github.com/project-fluxo/fluxo/blob/c7e0cc9b7fd4569dcab67bbb6e5a25c0a84859f1/src/equation/navierstokes/calctimestep.f90#L130-L132
             diffusivity1_transformed = diffusivity * (Ja11^2 + Ja12^2 + Ja13^2)
             Ja21, Ja22, Ja23 = get_contravariant_vector(2, contravariant_vectors,
                                                         i, j, k, element)
@@ -219,6 +225,12 @@ function max_dt(u, t,
             # Local diffusivity transformed to the reference element
             Ja11, Ja12, Ja13 = get_contravariant_vector(1, contravariant_vectors,
                                                         i, j, k, element)
+            # The metric terms are squared due to the second derivatives in the parabolic terms,
+            # which lead to application of the chain rule (coordinate transformation) twice.
+            # See for instance also the implementation in FLEXI
+            # https://github.com/flexi-framework/flexi/blob/e980e8635e6605daf906fc176c9897314a9cd9b1/src/equations/navierstokes/calctimestep.f90#L75-L77
+            # or FLUXO:
+            # https://github.com/project-fluxo/fluxo/blob/c7e0cc9b7fd4569dcab67bbb6e5a25c0a84859f1/src/equation/navierstokes/calctimestep.f90#L130-L132
             diffusivity1_transformed = diffusivity * (Ja11^2 + Ja12^2 + Ja13^2)
             Ja21, Ja22, Ja23 = get_contravariant_vector(2, contravariant_vectors,
                                                         i, j, k, element)

From 88909587026398173d7a94642538202c7830d5fe Mon Sep 17 00:00:00 2001
From: Daniel_Doehring <daniel.doehring@rwth-aachen.de>
Date: Wed, 10 Dec 2025 11:28:49 +0100
Subject: [PATCH 10/11] shorten test time

---
 test/test_parabolic_3d.jl | 23 +++++++++++------------
 1 file changed, 11 insertions(+), 12 deletions(-)

diff --git a/test/test_parabolic_3d.jl b/test/test_parabolic_3d.jl
index c5cc0885e18..1fbc038ead3 100644
--- a/test/test_parabolic_3d.jl
+++ b/test/test_parabolic_3d.jl
@@ -372,8 +372,7 @@ end
 @trixi_testset "P4estMesh3D: elixir_navierstokes_taylor_green_vortex.jl (Diffusive CFL)" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "p4est_3d_dgsem",
                                  "elixir_navierstokes_taylor_green_vortex.jl"),
-                        initial_refinement_level=2, tspan=(0.0, 0.25),
-                        surface_flux=FluxHLL(min_max_speed_naive),
+                        tspan=(0.0, 0.1),
                         mu=0.5, # render flow diffusion-dominated
                         callbacks=CallbackSet(summary_callback, analysis_callback,
                                               alive_callback,
@@ -382,18 +381,18 @@ end
                         adaptive=false, # respect CFL
                         ode_alg=CKLLSRK95_4S(),
                         l2=[
-                            0.0001336113719363628,
-                            0.1110290325891525,
-                            0.11102903258915098,
-                            0.009310736674102591,
-                            0.18713598666809547
+                            0.0001022410497625877,
+                            0.04954975879887512,
+                            0.049549758798875056,
+                            0.005853983721675305,
+                            0.09161121143324424
                         ],
                         linf=[
-                            0.0005497841741080034,
-                            0.31512261528156615,
-                            0.31512262804742686,
-                            0.02632785741783134,
-                            0.5709646222385913
+                            0.00039284994602417633,
+                            0.14026307274342587,
+                            0.14026307274350203,
+                            0.017003338595870714,
+                            0.2823457296549634
                         ])
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)

From 9752519c5443853c60f794bb104417fb155c48c4 Mon Sep 17 00:00:00 2001
From: Daniel_Doehring <daniel.doehring@rwth-aachen.de>
Date: Thu, 11 Dec 2025 15:02:41 +0100
Subject: [PATCH 11/11] try fixing linearsolve

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

diff --git a/test/Project.toml b/test/Project.toml
index 824aade8a6c..bb5b88b5649 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -55,7 +55,7 @@ FiniteDiff = "2.27.0"
 ForwardDiff = "0.10.36, 1"
 Krylov = "0.10"
 LinearAlgebra = "1"
-LinearSolve = "3.13"
+LinearSolve = "3.13, <3.49"
 MPI = "0.20.22"
 NLsolve = "4.5.1"
 OrdinaryDiffEqBDF = "1.1"