(Thread-)Parallelize bounds check routine for subcell IDP limiting

src/callbacks_stage/subcell_bounds_check.jl

@@ -118,7 +118,7 @@ end
 @inline function finalize_callback(callback::BoundsCheckCallback, semi,
                                    limiter::SubcellLimiterIDP)
     (; local_minmax, positivity) = limiter
-    (; idp_bounds_delta) = limiter.cache
+    (; idp_bounds_delta_threaded) = limiter.cache
     variables = varnames(cons2cons, semi.equations)
 
     println("─"^100)
@@ -128,8 +128,10 @@ end
         for v in limiter.local_minmax_variables_cons
             v_string = string(v)
             println("$(variables[v]):")
-            println("-lower bound: ", idp_bounds_delta[Symbol(v_string, "_min")][2])
-            println("-upper bound: ", idp_bounds_delta[Symbol(v_string, "_max")][2])
+            println("- lower bound: ",
+                    idp_bounds_delta_threaded[Symbol(v_string, "_min")][1][2])
+            println("- upper bound: ",
+                    idp_bounds_delta_threaded[Symbol(v_string, "_max")][1][2])
         end
     end
     if positivity
@@ -138,7 +140,7 @@ end
                 continue
             end
             println(string(variables[v]) * ":\n- positivity: ",
-                    idp_bounds_delta[Symbol(string(v), "_min")][2])
+                    idp_bounds_delta_threaded[Symbol(string(v), "_min")][1][2])
         end
     end
     println("─"^100 * "\n")

src/callbacks_stage/subcell_bounds_check_2d.jl

@@ -10,26 +10,38 @@
                               time, iter, output_directory, save_errors)
     (; local_minmax, positivity) = solver.volume_integral.limiter
     (; variable_bounds) = limiter.cache.subcell_limiter_coefficients
-    (; idp_bounds_delta) = limiter.cache
+    (; idp_bounds_delta_threaded) = limiter.cache
 
     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 = idp_bounds_delta[key_min]
-            deviation_max = idp_bounds_delta[key_max]
-            for element in eachelement(solver, cache), j in eachnode(solver),
-                i in eachnode(solver)
-
-                var = u[v, i, j, element]
-                deviation_min[1] = max(deviation_min[1],
-                                       variable_bounds[key_min][i, j, element] - var)
-                deviation_max[1] = max(deviation_max[1],
-                                       var - variable_bounds[key_max][i, j, element])
+            deviation_min_threaded = idp_bounds_delta_threaded[key_min]
+            deviation_max_threaded = idp_bounds_delta_threaded[key_max]
+            @threaded for element in eachelement(solver, cache)
+                deviation_min = deviation_min_threaded[Threads.threadid()]
+                deviation_max = deviation_max_threaded[Threads.threadid()]
+                for j in eachnode(solver), i in eachnode(solver)
+                    var = u[v, i, j, element]
+                    deviation_min[1] = max(deviation_min[1],
+                                           variable_bounds[key_min][i, j, element] -
+                                           var)
+                    deviation_max[1] = max(deviation_max[1],
+                                           var -
+                                           variable_bounds[key_max][i, j, element])
+                end
             end
-            deviation_min[2] = max(deviation_min[2], deviation_min[1])
-            deviation_max[2] = max(deviation_max[2], deviation_max[1])
+            for i in 2:Threads.nthreads()
+                deviation_min_threaded[1][1] = max(deviation_min_threaded[1][1],
+                                                   deviation_min_threaded[i][1])
+                deviation_max_threaded[1][1] = max(deviation_max_threaded[1][1],
+                                                   deviation_max_threaded[i][1])
+            end
+            deviation_min_threaded[1][2] = max(deviation_min_threaded[1][2],
+                                               deviation_min_threaded[1][1])
+            deviation_max_threaded[1][2] = max(deviation_max_threaded[1][2],
+                                               deviation_max_threaded[1][1])
         end
     end
     if positivity
@@ -38,15 +50,21 @@
                 continue
             end
             key = Symbol(string(v), "_min")
-            deviation = idp_bounds_delta[key]
-            for element in eachelement(solver, cache), j in eachnode(solver),
-                i in eachnode(solver)
-
-                var = u[v, i, j, element]
-                deviation[1] = max(deviation[1],
-                                   variable_bounds[key][i, j, element] - var)
+            deviation_threaded = idp_bounds_delta_threaded[key]
+            @threaded for element in eachelement(solver, cache)
+                deviation = deviation_threaded[Threads.threadid()]
+                for j in eachnode(solver), i in eachnode(solver)
+                    var = u[v, i, j, element]
+                    deviation[1] = max(deviation[1],
+                                       variable_bounds[key][i, j, element] - var)
+                end
+            end
+            for i in 2:Threads.nthreads()
+                deviation_threaded[1][1] = max(deviation_threaded[1][1],
+                                               deviation_threaded[i][1])
             end
-            deviation[2] = max(deviation[2], deviation[1])
+            deviation_threaded[1][2] = max(deviation_threaded[1][2],
+                                           deviation_threaded[1][1])
         end
     end
     if save_errors
@@ -56,23 +74,28 @@
             if local_minmax
                 for v in limiter.local_minmax_variables_cons
                     v_string = string(v)
-                    print(f, ", ", idp_bounds_delta[Symbol(v_string, "_min")][1], ", ",
-                          idp_bounds_delta[Symbol(v_string, "_max")][1])
+                    print(f, ", ",
+                          idp_bounds_delta_threaded[Symbol(v_string, "_min")][1][1],
+                          ", ",
+                          idp_bounds_delta_threaded[Symbol(v_string, "_max")][1][1])
                 end
             end
             if positivity
                 for v in limiter.positivity_variables_cons
                     if v in limiter.local_minmax_variables_cons
                         continue
                     end
-                    print(f, ", ", idp_bounds_delta[Symbol(string(v), "_min")][1])
+                    print(f, ", ",
+                          idp_bounds_delta_threaded[Symbol(string(v), "_min")][1][1])
                 end
             end
             println(f)
         end
-        # Reset first entries of idp_bounds_delta
-        for (key, _) in idp_bounds_delta
-            idp_bounds_delta[key][1] = zero(eltype(idp_bounds_delta[key][1]))
+        # Reset first entries of idp_bounds_delta_threaded
+        for (key, _) in idp_bounds_delta_threaded
+            for i in 1:Threads.nthreads()
+                idp_bounds_delta_threaded[key][i][1] = zero(eltype(idp_bounds_delta_threaded[key][i][1]))
+            end
         end
     end
 

src/solvers/dgsem_tree/subcell_limiters_2d.jl

@@ -13,22 +13,23 @@ function create_cache(limiter::Type{SubcellLimiterIDP}, equations::AbstractEquat
                                                                         nnodes(basis),
                                                                         bound_keys)
 
-    # Memory for bounds checking routine with `BoundsCheckCallback`.
+    # Threaded memory for bounds checking routine with `BoundsCheckCallback`.
     # The first entry of each vector contains the maximum deviation since the last export.
-    # The second one contains the total maximum deviation.
-    idp_bounds_delta = Dict{Symbol, Vector{real(basis)}}()
+    # In the second entry, the total maximum deviation is saved.
+    idp_bounds_delta_threaded = Dict{Symbol, Vector{Vector{real(basis)}}}()
     for key in bound_keys
-        idp_bounds_delta[key] = zeros(real(basis), 2)
+        idp_bounds_delta_threaded[key] = [zeros(real(basis), 2)
+                                          for _ in 1:Threads.nthreads()]
     end
 
-    return (; subcell_limiter_coefficients, idp_bounds_delta)
+    return (; subcell_limiter_coefficients, idp_bounds_delta_threaded)
 end
 
 function (limiter::SubcellLimiterIDP)(u::AbstractArray{<:Any, 4}, semi, dg::DGSEM, t,
                                       dt;
                                       kwargs...)
     @unpack alpha = limiter.cache.subcell_limiter_coefficients
-    alpha .= zero(eltype(alpha))
+    @trixi_timeit timer() "reset alpha" reset_du!(alpha, dg, semi.cache)
 
     if limiter.local_minmax
         @trixi_timeit timer() "local min/max limiting" idp_local_minmax!(alpha, limiter,