Merge pull request #117 from bennibolm/subcell-limiting-bounds-check-fix

Revise calculation of bound check; Replace construction of Symbols

src/callbacks_stage/subcell_bounds_check.jl

@@ -86,7 +86,7 @@ function init_callback(callback::BoundsCheckCallback, semi, limiter::SubcellLimi
         print(f, "# iter, simu_time")
         if positivity
             for v in limiter.positivity_variables_cons
-                print(f, ", $(variables[v])_min")
+                print(f, ", " * string(variables[v]) * "_min")
             end
         end
         println(f)
@@ -117,8 +117,8 @@ end
     println("─"^100)
     if positivity
         for v in limiter.positivity_variables_cons
-            println("$(variables[v]):\n- positivity: ",
-                    idp_bounds_delta[Symbol("$(v)_min")])
+            println(string(variables[v]) * ":\n- positivity: ",
+                    idp_bounds_delta[Symbol(string(v), "_min")][2])
         end
     end
     println("─"^100 * "\n")

src/callbacks_stage/subcell_bounds_check_2d.jl

@@ -12,35 +12,36 @@
     (; variable_bounds) = limiter.cache.subcell_limiter_coefficients
     (; idp_bounds_delta) = limiter.cache
 
-    if save_errors
-        open("$output_directory/deviations.txt", "a") do f
-            print(f, iter, ", ", time)
-        end
-    end
     if positivity
         for v in limiter.positivity_variables_cons
-            key = Symbol("$(v)_min")
-            deviation_min = zero(eltype(u))
+            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_min = max(deviation_min,
-                                    variable_bounds[key][i, j, element] - var)
-            end
-            idp_bounds_delta[key] = max(idp_bounds_delta[key], deviation_min)
-            if save_errors
-                deviation_min_ = deviation_min
-                open("$output_directory/deviations.txt", "a") do f
-                    print(f, ", ", deviation_min_)
-                end
+                deviation[1] = max(deviation[1],
+                                   variable_bounds[key][i, j, element] - var)
             end
+            deviation[2] = max(deviation[2], deviation[1])
         end
     end
     if save_errors
+        # Print to output file
         open("$output_directory/deviations.txt", "a") do f
+            print(f, iter, ", ", time)
+            if positivity
+                for v in limiter.positivity_variables_cons
+                    key = Symbol(string(v), "_min")
+                    print(f, ", ", idp_bounds_delta[key][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]))
+        end
     end
 
     return nothing

src/solvers/dgsem_tree/subcell_limiters.jl

@@ -53,7 +53,7 @@ function SubcellLimiterIDP(equations::AbstractEquations, basis;
                            positivity_correction_factor = 0.1)
     positivity = (length(positivity_variables_cons) > 0)
 
-    bound_keys = Tuple(Symbol("$(i)_min") for i in positivity_variables_cons)
+    bound_keys = Tuple(Symbol(string(v), "_min") for v in positivity_variables_cons)
 
     cache = create_cache(SubcellLimiterIDP, equations, basis, bound_keys)
 

src/solvers/dgsem_tree/subcell_limiters_2d.jl

@@ -13,9 +13,12 @@ function create_cache(limiter::Type{SubcellLimiterIDP}, equations::AbstractEquat
                                                                         nnodes(basis),
                                                                         bound_keys)
 
-    idp_bounds_delta = Dict{Symbol, real(basis)}()
+    # 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)}}()
     for key in bound_keys
-        idp_bounds_delta[key] = zero(real(basis))
+        idp_bounds_delta[key] = zeros(real(basis), 2)
     end
 
     return (; subcell_limiter_coefficients, idp_bounds_delta)
@@ -65,7 +68,7 @@ end
     (; positivity_correction_factor) = limiter
 
     (; variable_bounds) = limiter.cache.subcell_limiter_coefficients
-    var_min = variable_bounds[Symbol("$(variable)_min")]
+    var_min = variable_bounds[Symbol(string(variable), "_min")]
 
     @threaded for element in eachelement(dg, semi.cache)
         inverse_jacobian = cache.elements.inverse_jacobian[element]