diff --git a/src/solvers/dgsem_tree/subcell_limiters.jl b/src/solvers/dgsem_tree/subcell_limiters.jl
index 554c4810178..86d70e3853a 100644
--- a/src/solvers/dgsem_tree/subcell_limiters.jl
+++ b/src/solvers/dgsem_tree/subcell_limiters.jl
@@ -25,16 +25,16 @@ end
 Subcell invariant domain preserving (IDP) limiting used with [`VolumeIntegralSubcellLimiting`](@ref)
 including:
 - Local maximum/minimum Zalesak-type limiting for conservative variables (`local_minmax_variables_cons`)
-- Positivity limiting for conservative variables (`positivity_variables_cons`) and non-linear variables
+- Positivity limiting for conservative variables (`positivity_variables_cons`) and nonlinear variables
 (`positivity_variables_nonlinear`)
 
 Conservative variables to be limited are passed as a vector of strings, e.g. `local_minmax_variables_cons = ["rho"]`
-and `positivity_variables_cons = ["rho"]`. For non-linear variables the specific functions are
+and `positivity_variables_cons = ["rho"]`. For nonlinear variables the specific functions are
 passed in a vector, e.g. `positivity_variables_nonlinear = [pressure]`.
 
 The bounds are calculated using the low-order FV solution. The positivity limiter uses
 `positivity_correction_factor` such that `u^new >= positivity_correction_factor * u^FV`.
-The limiting of non-linear variables uses a Newton-bisection method with a maximum of
+The limiting of nonlinear variables uses a Newton-bisection method with a maximum of
 `max_iterations_newton` iterations, relative and absolute tolerances of `newton_tolerances`
 and a gamma constant of `gamma_constant_newton` (`gamma_constant_newton>=2*d`,
 where `d = #dimensions`).
@@ -124,9 +124,9 @@ function Base.show(io::IO, limiter::SubcellLimiterIDP)
     if !(local_minmax || positivity)
         print(io, "No limiter selected => pure DG method")
     else
-        print(io, "limiter=(")
-        local_minmax && print(io, "min/max limiting, ")
-        positivity && print(io, "positivity")
+        print(io, "Limiter=(")
+        local_minmax && print(io, "Local min/max, ")
+        positivity && print(io, "Positivity, ")
         print(io, "), ")
     end
     print(io, "Local bounds with FV solution")
@@ -147,15 +147,15 @@ function Base.show(io::IO, ::MIME"text/plain", limiter::SubcellLimiterIDP)
             if local_minmax
                 setup = [
                     setup...,
-                    "" => "local maximum/minimum bounds for conservative variables $(limiter.local_minmax_variables_cons)",
+                    "" => "Local maximum/minimum limiting for conservative variables $(limiter.local_minmax_variables_cons)",
                 ]
             end
             if positivity
-                string = "positivity for conservative variables $(limiter.positivity_variables_cons) and $(limiter.positivity_variables_nonlinear)"
+                string = "Positivity limiting for conservative variables $(limiter.positivity_variables_cons) and $(limiter.positivity_variables_nonlinear)"
                 setup = [setup..., "" => string]
                 setup = [
                     setup...,
-                    "" => "   positivity correction factor = $(limiter.positivity_correction_factor)",
+                    "" => "- with positivity correction factor = $(limiter.positivity_correction_factor)",
                 ]
             end
             setup = [
diff --git a/src/solvers/dgsem_tree/subcell_limiters_2d.jl b/src/solvers/dgsem_tree/subcell_limiters_2d.jl
index e362d71cacc..1b6e8d318e7 100644
--- a/src/solvers/dgsem_tree/subcell_limiters_2d.jl
+++ b/src/solvers/dgsem_tree/subcell_limiters_2d.jl
@@ -5,6 +5,10 @@
 @muladd begin
 #! format: noindent
 
+###############################################################################
+# IDP Limiting
+###############################################################################
+
 # this method is used when the limiter is constructed as for shock-capturing volume integrals
 function create_cache(limiter::Type{SubcellLimiterIDP}, equations::AbstractEquations{2},
                       basis::LobattoLegendreBasis, bound_keys)
@@ -56,6 +60,9 @@ function (limiter::SubcellLimiterIDP)(u::AbstractArray{<:Any, 4}, semi, dg::DGSE
     return nothing
 end
 
+###############################################################################
+# Calculation of local bounds using low-order FV solution
+
 @inline function calc_bounds_twosided!(var_min, var_max, variable, u, t, semi)
     mesh, equations, dg, cache = mesh_equations_solver_cache(semi)
     # Calc bounds inside elements
@@ -154,6 +161,9 @@ end
     return nothing
 end
 
+###############################################################################
+# Local minimum/maximum limiting
+
 @inline function idp_local_minmax!(alpha, limiter, u, t, dt, semi)
     for variable in limiter.local_minmax_variables_cons
         idp_local_minmax!(alpha, limiter, u, t, dt, semi, variable)
@@ -223,20 +233,32 @@ end
     return nothing
 end
 
+###############################################################################
+# Global positivity limiting
+
 @inline function idp_positivity!(alpha, limiter, u, dt, semi)
     # Conservative variables
     for variable in limiter.positivity_variables_cons
-        idp_positivity!(alpha, limiter, u, dt, semi, variable)
+        @trixi_timeit timer() "conservative variables" idp_positivity!(alpha, limiter,
+                                                                       u, dt, semi,
+                                                                       variable)
     end
 
     # Nonlinear variables
     for variable in limiter.positivity_variables_nonlinear
-        idp_positivity_nonlinear!(alpha, limiter, u, dt, semi, variable)
+        @trixi_timeit timer() "nonlinear variables" idp_positivity_nonlinear!(alpha,
+                                                                              limiter,
+                                                                              u, dt,
+                                                                              semi,
+                                                                              variable)
     end
 
     return nothing
 end
 
+###############################################################################
+# Global positivity limiting of conservative variables
+
 @inline function idp_positivity!(alpha, limiter, u, dt, semi, variable)
     mesh, equations, dg, cache = mesh_equations_solver_cache(semi)
     (; antidiffusive_flux1_L, antidiffusive_flux2_L, antidiffusive_flux1_R, antidiffusive_flux2_R) = cache.antidiffusive_fluxes
@@ -299,13 +321,14 @@ end
 end
 
 @inline function idp_positivity_nonlinear!(alpha, limiter, u, dt, semi, variable)
-    mesh, equations, dg, cache = mesh_equations_solver_cache(semi)
+    _, equations, dg, cache = mesh_equations_solver_cache(semi)
     (; positivity_correction_factor) = limiter
 
     (; variable_bounds) = limiter.cache.subcell_limiter_coefficients
     var_min = variable_bounds[Symbol(string(variable), "_min")]
 
     @threaded for element in eachelement(dg, semi.cache)
+        inverse_jacobian = cache.elements.inverse_jacobian[element]
         for j in eachnode(dg), i in eachnode(dg)
             # Compute bound
             u_local = get_node_vars(u, equations, dg, i, j, element)
@@ -318,8 +341,8 @@ end
             # Perform Newton's bisection method to find new alpha
             newton_loops_alpha!(alpha, var_min[i, j, element], u_local, i, j, element,
                                 variable, initial_check_nonnegative,
-                                final_check_nonnegative,
-                                dt, mesh, equations, dg, cache, limiter)
+                                final_check_nonnegative, inverse_jacobian,
+                                dt, equations, dg, cache, limiter)
         end
     end
 
@@ -327,11 +350,10 @@ end
 end
 
 @inline function newton_loops_alpha!(alpha, bound, u, i, j, element, variable,
-                                     initial_check, final_check, dt, mesh, equations,
-                                     dg, cache, limiter)
+                                     initial_check, final_check, inverse_jacobian, dt,
+                                     equations, dg, cache, limiter)
     (; inverse_weights) = dg.basis
     (; antidiffusive_flux1_L, antidiffusive_flux2_L, antidiffusive_flux1_R, antidiffusive_flux2_R) = cache.antidiffusive_fluxes
-    inverse_jacobian = cache.elements.inverse_jacobian[element]
 
     (; gamma_constant_newton) = limiter