Revise parallel memory structure

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_threaded) = limiter.cache
+    (; idp_bounds_delta_global) = limiter.cache
     variables = varnames(cons2cons, semi.equations)
 
     println("─"^100)
@@ -129,9 +129,9 @@ end
             v_string = string(v)
             println("$(variables[v]):")
             println("- lower bound: ",
-                    idp_bounds_delta_threaded[Symbol(v_string, "_min")][1][2])
+                    idp_bounds_delta_global[Symbol(v_string, "_min")])
             println("- upper bound: ",
-                    idp_bounds_delta_threaded[Symbol(v_string, "_max")][1][2])
+                    idp_bounds_delta_global[Symbol(v_string, "_max")])
         end
     end
     if positivity
@@ -140,7 +140,7 @@ end
                 continue
             end
             println(string(variables[v]) * ":\n- positivity: ",
-                    idp_bounds_delta_threaded[Symbol(string(v), "_min")][1][2])
+                    idp_bounds_delta_global[Symbol(string(v), "_min")])
         end
     end
     println("─"^100 * "\n")

src/callbacks_stage/subcell_bounds_check_2d.jl

@@ -10,27 +10,27 @@
                               time, iter, output_directory, save_errors)
     (; local_minmax, positivity) = solver.volume_integral.limiter
     (; variable_bounds) = limiter.cache.subcell_limiter_coefficients
-    (; idp_bounds_delta_threaded) = limiter.cache
+    (; idp_bounds_delta_local, idp_bounds_delta_global) = 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_threaded = idp_bounds_delta_threaded[key_min]
-            deviation_max_threaded = idp_bounds_delta_threaded[key_max]
+            deviation_min_threaded = idp_bounds_delta_local[key_min]
+            deviation_max_threaded = idp_bounds_delta_local[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])
+                    deviation_min = max(deviation_min,
+                                        variable_bounds[key_min][i, j, element] - var)
+                    deviation_max = max(deviation_max,
+                                        var - variable_bounds[key_max][i, j, element])
                 end
+                deviation_min_threaded[Threads.threadid()] = deviation_min
+                deviation_max_threaded[Threads.threadid()] = deviation_max
             end
         end
     end
@@ -40,27 +40,25 @@
                 continue
             end
             key = Symbol(string(v), "_min")
-            deviation_threaded = idp_bounds_delta_threaded[key]
+            deviation_threaded = idp_bounds_delta_local[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)
+                    deviation = max(deviation,
+                                    variable_bounds[key][i, j, element] - var)
                 end
+                deviation_threaded[Threads.threadid()] = deviation
             end
         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
+    for (key, _) in idp_bounds_delta_local
+        # Reduce threaded local maximum deviations. Save in first entry.
+        idp_bounds_delta_local[key][1] = reduce(max, idp_bounds_delta_local[key])
         # 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])
+        idp_bounds_delta_global[key] = max(idp_bounds_delta_global[key],
+                                           idp_bounds_delta_local[key][1])
     end
 
     if save_errors
@@ -71,27 +69,23 @@
                 for v in limiter.local_minmax_variables_cons
                     v_string = string(v)
                     print(f, ", ",
-                          idp_bounds_delta_threaded[Symbol(v_string, "_min")][1][1],
-                          ", ",
-                          idp_bounds_delta_threaded[Symbol(v_string, "_max")][1][1])
+                          idp_bounds_delta_local[Symbol(v_string, "_min")][1], ", ",
+                          idp_bounds_delta_local[Symbol(v_string, "_max")][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_threaded[Symbol(string(v), "_min")][1][1])
+                    print(f, ", ", idp_bounds_delta_local[Symbol(string(v), "_min")][1])
                 end
             end
             println(f)
         end
-        # 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
+        # Reset local maximum deviations
+        for (key, _) in idp_bounds_delta_local
+            idp_bounds_delta_local[key] .= zero(eltype(idp_bounds_delta_local[key][1]))
         end
     end
 

src/solvers/dgsem_tree/subcell_limiters_2d.jl

@@ -13,16 +13,19 @@ function create_cache(limiter::Type{SubcellLimiterIDP}, equations::AbstractEquat
                                                                         nnodes(basis),
                                                                         bound_keys)
 
-    # Threaded memory for bounds checking routine with `BoundsCheckCallback`.
-    # The first entry of each vector contains the maximum deviation since the last export.
-    # In the second entry, the total maximum deviation is saved.
-    idp_bounds_delta_threaded = Dict{Symbol, Vector{Vector{real(basis)}}}()
+    # Memory for bounds checking routine with `BoundsCheckCallback`.
+    # Local variable contains the maximum deviation since the last export.
+    # Using threaded variable to parallelize bounds check.
+    idp_bounds_delta_local = Dict{Symbol, Vector{real(basis)}}()
+    # Global variable contains the total maximum deviation.
+    idp_bounds_delta_global = Dict{Symbol, real(basis)}()
     for key in bound_keys
-        idp_bounds_delta_threaded[key] = [zeros(real(basis), 2)
-                                          for _ in 1:Threads.nthreads()]
+        idp_bounds_delta_local[key] = [zero(real(basis)) for _ in 1:Threads.nthreads()]
+        idp_bounds_delta_global[key] = zero(real(basis))
     end
 
-    return (; subcell_limiter_coefficients, idp_bounds_delta_threaded)
+    return (; subcell_limiter_coefficients, idp_bounds_delta_local,
+            idp_bounds_delta_global)
 end
 
 function (limiter::SubcellLimiterIDP)(u::AbstractArray{<:Any, 4}, semi, dg::DGSEM, t,