Add bounds check for local limiting

examples/tree_2d_dgsem/elixir_euler_blast_wave_sc_subcell.jl

@@ -88,7 +88,7 @@ callbacks = CallbackSet(summary_callback,
 ###############################################################################
 # run the simulation
 
-stage_callbacks = (SubcellLimiterIDPCorrection(),)
+stage_callbacks = (SubcellLimiterIDPCorrection(), BoundsCheckCallback(save_errors=false))
 
 sol = Trixi.solve(ode, Trixi.SimpleSSPRK33(stage_callbacks=stage_callbacks);
                   dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback

examples/tree_2d_dgsem/elixir_euler_sedov_blast_wave_sc_subcell.jl

@@ -84,7 +84,7 @@ callbacks = CallbackSet(summary_callback,
 ###############################################################################
 # run the simulation
 
-stage_callbacks = (SubcellLimiterIDPCorrection(),)
+stage_callbacks = (SubcellLimiterIDPCorrection(), BoundsCheckCallback(save_errors=false))
 
 sol = Trixi.solve(ode, Trixi.SimpleSSPRK33(stage_callbacks=stage_callbacks);
                   dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback

examples/tree_2d_dgsem/elixir_eulermulti_shock_bubble_shockcapturing_subcell_minmax.jl

@@ -131,7 +131,7 @@ callbacks           = CallbackSet(summary_callback,
 ###############################################################################
 # run the simulation
 
-stage_callbacks = (SubcellLimiterIDPCorrection(),)
+stage_callbacks = (SubcellLimiterIDPCorrection(), BoundsCheckCallback(save_errors=false))
 
 sol = Trixi.solve(ode, Trixi.SimpleSSPRK33(stage_callbacks=stage_callbacks);
                   dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback

src/callbacks_stage/subcell_bounds_check.jl

@@ -77,15 +77,24 @@ function init_callback(callback::BoundsCheckCallback, semi, limiter::SubcellLimi
         return nothing
     end
 
-    (; positivity) = limiter
+    (; local_minmax, positivity) = limiter
     (; output_directory) = callback
     variables = varnames(cons2cons, semi.equations)
 
     mkpath(output_directory)
     open("$output_directory/deviations.txt", "a") do f
         print(f, "# iter, simu_time")
+        if local_minmax
+            for v in limiter.local_minmax_variables_cons
+                variable_string = string(variables[v])
+                print(f, ", " * variable_string * "_min, " * variable_string * "_max")
+            end
+        end
         if positivity
             for v in limiter.positivity_variables_cons
+                if v in limiter.local_minmax_variables_cons
+                    continue
+                end
                 print(f, ", " * string(variables[v]) * "_min")
             end
         end
@@ -108,15 +117,26 @@ end
 
 @inline function finalize_callback(callback::BoundsCheckCallback, semi,
                                    limiter::SubcellLimiterIDP)
-    (; positivity) = limiter
+    (; local_minmax, positivity) = limiter
     (; idp_bounds_delta) = limiter.cache
     variables = varnames(cons2cons, semi.equations)
 
     println("─"^100)
     println("Maximum deviation from bounds:")
     println("─"^100)
+    if local_minmax
+        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])
+        end
+    end
     if positivity
         for v in limiter.positivity_variables_cons
+            if v in limiter.local_minmax_variables_cons
+                continue
+            end
             println(string(variables[v]) * ":\n- positivity: ",
                     idp_bounds_delta[Symbol(string(v), "_min")][2])
         end

src/callbacks_stage/subcell_bounds_check_2d.jl

@@ -8,12 +8,35 @@
 @inline function check_bounds(u, mesh::AbstractMesh{2}, equations, solver, cache,
                               limiter::SubcellLimiterIDP,
                               time, iter, output_directory, save_errors)
-    (; positivity) = solver.volume_integral.limiter
+    (; local_minmax, positivity) = solver.volume_integral.limiter
     (; variable_bounds) = limiter.cache.subcell_limiter_coefficients
     (; idp_bounds_delta) = 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])
+            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
         for v in limiter.positivity_variables_cons
+            if v in limiter.local_minmax_variables_cons
+                continue
+            end
             key = Symbol(string(v), "_min")
             deviation = idp_bounds_delta[key]
             for element in eachelement(solver, cache), j in eachnode(solver),
@@ -30,10 +53,19 @@
         # 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])
+                end
+            end
             if positivity
                 for v in limiter.positivity_variables_cons
-                    key = Symbol(string(v), "_min")
-                    print(f, ", ", idp_bounds_delta[key][1])
+                    if v in limiter.local_minmax_variables_cons
+                        continue
+                    end
+                    print(f, ", ", idp_bounds_delta[Symbol(string(v), "_min")][1])
                 end
             end
             println(f)

src/solvers/dgsem_tree/subcell_limiters.jl

@@ -62,8 +62,9 @@ function SubcellLimiterIDP(equations::AbstractEquations, basis;
     bound_keys = ()
     if local_minmax
         for v in local_minmax_variables_cons
-            bound_keys = (bound_keys..., Symbol(string(v), "_min"),
-                          Symbol(string(v), "_max"))
+            v_string = string(v)
+            bound_keys = (bound_keys..., Symbol(v_string, "_min"),
+                          Symbol(v_string, "_max"))
         end
     end
     for v in positivity_variables_cons