diff --git a/src/solvers/dgsem_tree/subcell_limiters_2d.jl b/src/solvers/dgsem_tree/subcell_limiters_2d.jl
index 29a1ac322fa..832140e4ff4 100644
--- a/src/solvers/dgsem_tree/subcell_limiters_2d.jl
+++ b/src/solvers/dgsem_tree/subcell_limiters_2d.jl
@@ -17,30 +17,33 @@
     mesh, _, dg, cache = mesh_equations_solver_cache(semi)
     # Calc bounds inside elements
     @threaded for element in eachelement(dg, cache)
-        var_min[:, :, element] .= typemax(eltype(var_min))
-        var_max[:, :, element] .= typemin(eltype(var_max))
-        # Calculate bounds at Gauss-Lobatto nodes using u
+        # Calculate bounds at Gauss-Lobatto nodes
         for j in eachnode(dg), i in eachnode(dg)
             var = u[variable, i, j, element]
+            var_min[i, j, element] = var
+            var_max[i, j, element] = var
+        end
+
+        # Apply values in x direction
+        for j in eachnode(dg), i in 2:nnodes(dg)
+            var = u[variable, i - 1, j, element]
             var_min[i, j, element] = min(var_min[i, j, element], var)
             var_max[i, j, element] = max(var_max[i, j, element], var)
 
-            if i > 1
-                var_min[i - 1, j, element] = min(var_min[i - 1, j, element], var)
-                var_max[i - 1, j, element] = max(var_max[i - 1, j, element], var)
-            end
-            if i < nnodes(dg)
-                var_min[i + 1, j, element] = min(var_min[i + 1, j, element], var)
-                var_max[i + 1, j, element] = max(var_max[i + 1, j, element], var)
-            end
-            if j > 1
-                var_min[i, j - 1, element] = min(var_min[i, j - 1, element], var)
-                var_max[i, j - 1, element] = max(var_max[i, j - 1, element], var)
-            end
-            if j < nnodes(dg)
-                var_min[i, j + 1, element] = min(var_min[i, j + 1, element], var)
-                var_max[i, j + 1, element] = max(var_max[i, j + 1, element], var)
-            end
+            var = u[variable, i, j, element]
+            var_min[i - 1, j, element] = min(var_min[i - 1, j, element], var)
+            var_max[i - 1, j, element] = max(var_max[i - 1, j, element], var)
+        end
+
+        # Apply values in y direction
+        for j in 2:nnodes(dg), i in eachnode(dg)
+            var = u[variable, i, j - 1, element]
+            var_min[i, j, element] = min(var_min[i, j, element], var)
+            var_max[i, j, element] = max(var_max[i, j, element], var)
+
+            var = u[variable, i, j, element]
+            var_min[i, j - 1, element] = min(var_min[i, j - 1, element], var)
+            var_max[i, j - 1, element] = max(var_max[i, j - 1, element], var)
         end
     end
 
@@ -120,6 +123,10 @@ end
 @inline function calc_bounds_onesided!(var_minmax, min_or_max, variable,
                                        u::AbstractArray{<:Any, 4}, t, semi)
     mesh, equations, dg, cache = mesh_equations_solver_cache(semi)
+
+    # The approach used in `calc_bounds_twosided!` is not used here because it requires more
+    # evaluations of the variable and is therefore slower.
+
     # Calc bounds inside elements
     @threaded for element in eachelement(dg, cache)
         # Reset bounds
@@ -131,7 +138,7 @@ end
             end
         end
 
-        # Calculate bounds at Gauss-Lobatto nodes using u
+        # Calculate bounds at Gauss-Lobatto nodes
         for j in eachnode(dg), i in eachnode(dg)
             var = variable(get_node_vars(u, equations, dg, i, j, element), equations)
             var_minmax[i, j, element] = min_or_max(var_minmax[i, j, element], var)