(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,28 @@
                               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])
         end
     end
     if positivity
@@ -38,41 +40,58 @@
                 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
-            deviation[2] = max(deviation[2], deviation[1])
         end
     end
+
+    for (key, _) in idp_bounds_delta_threaded
+        # Calculate maximum deviations of all threads
+        for i in 2:Threads.nthreads()
+            idp_bounds_delta_threaded[key][1][1] = max(idp_bounds_delta_threaded[key][1][1],
+                                                       idp_bounds_delta_threaded[key][i][1])
+        end
+        # Update global maximum deviations
+        idp_bounds_delta_threaded[key][1][2] = max(idp_bounds_delta_threaded[key][1][2],
+                                                   idp_bounds_delta_threaded[key][1][1])
+    end
+
     if save_errors
         # Print to output file
         open("$output_directory/deviations.txt", "a") do f
             print(f, iter, ", ", time)
             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/dg_2d_subcell_limiters.jl

@@ -470,6 +470,9 @@ end
 For subcell limiting, the calculation of local bounds for non-periodic domains require the boundary
 outer state. This function returns the boundary value at time `t` and for node with spatial
 indices `indices`.
+
+!!! warning "Experimental implementation"
+    This is an experimental feature and may change in future releases.
 """
 @inline function get_boundary_outer_state(boundary_condition::BoundaryConditionDirichlet,
                                           cache, t, equations, dg, indices...)

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,