Add implementation of subcell bounds check

src/callbacks_stage/subcell_bounds_check.jl

@@ -8,9 +8,14 @@
 """
     BoundsCheckCallback(; output_directory="out", save_errors=false, interval=1)
 
-Bounds checking routine for [`SubcellLimiterIDP`](@ref). Applied as a stage callback for SSPRK
-methods. If `save_errors` is `true`, the resulting deviations are saved in
-`output_directory/deviations.txt` for every `interval` time steps.
+Subcell limiting techniques with [`SubcellLimiterIDP`](@ref) are constructed to adhere certain
+local or global bounds. To make sure that these bounds are actually met, this callback calculates
+the maximum deviation from the bounds. The maximum deviation per applied bound is printed to
+the screen at the end of the simulation.
+Additionally, for more insights, the occurring errors can be exported with `save_errors==true`
+every `interval` time steps during the simulation. Then, the maximum deviations since the last
+export are saved in "`output_directory`/deviations.txt".
+It has to be implied as a stage callback for SSPRK.
 
 !!! note
     For `SubcellLimiterIDP`, the solution is corrected in the a posteriori correction stage
@@ -33,14 +38,14 @@ function (callback::BoundsCheckCallback)(u_ode, integrator, stage)
     (; t, iter, alg) = integrator
     u = wrap_array(u_ode, mesh, equations, solver, cache)
 
-    save_errors_ = callback.save_errors && (callback.interval > 0) &&
-                   (stage == length(alg.c)) &&
-                   (iter % callback.interval == 0 || integrator.finalstep)
+    save_errors = callback.save_errors && (callback.interval > 0) &&
+                  (stage == length(alg.c)) &&
+                  (iter % callback.interval == 0 || integrator.finalstep)
     @trixi_timeit timer() "check_bounds" check_bounds(u, mesh, equations, solver, cache,
                                                       solver.volume_integral, t,
                                                       iter + 1,
                                                       callback.output_directory,
-                                                      save_errors_)
+                                                      save_errors)
 end
 
 function check_bounds(u, mesh, equations, solver, cache,
@@ -56,13 +61,14 @@ function check_bounds(u, mesh, equations, solver, cache,
                  output_directory, save_errors)
 end
 
-function init_callback(callback, semi)
+function init_callback(callback::BoundsCheckCallback, semi)
     init_callback(callback, semi, semi.solver.volume_integral)
 end
 
-init_callback(callback, semi, volume_integral::AbstractVolumeIntegral) = nothing
+init_callback(callback::BoundsCheckCallback, semi, volume_integral::AbstractVolumeIntegral) = nothing
 
-function init_callback(callback, semi, volume_integral::VolumeIntegralSubcellLimiting)
+function init_callback(callback::BoundsCheckCallback, semi,
+                       volume_integral::VolumeIntegralSubcellLimiting)
     init_callback(callback, semi, volume_integral.limiter)
 end
 
@@ -79,8 +85,8 @@ function init_callback(callback::BoundsCheckCallback, semi, limiter::SubcellLimi
     open("$output_directory/deviations.txt", "a") do f
         print(f, "# iter, simu_time")
         if positivity
-            for index in limiter.positivity_variables_cons
-                print(f, ", $(variables[index])_min")
+            for v in limiter.positivity_variables_cons
+                print(f, ", " * string(variables[v]) * "_min")
             end
         end
         println(f)
@@ -89,13 +95,13 @@ function init_callback(callback::BoundsCheckCallback, semi, limiter::SubcellLimi
     return nothing
 end
 
-function finalize_callback(callback, semi)
+function finalize_callback(callback::BoundsCheckCallback, semi)
     finalize_callback(callback, semi, semi.solver.volume_integral)
 end
 
-finalize_callback(callback, semi, volume_integral::AbstractVolumeIntegral) = nothing
+finalize_callback(callback::BoundsCheckCallback, semi, volume_integral::AbstractVolumeIntegral) = nothing
 
-function finalize_callback(callback, semi,
+function finalize_callback(callback::BoundsCheckCallback, semi,
                            volume_integral::VolumeIntegralSubcellLimiting)
     finalize_callback(callback, semi, volume_integral.limiter)
 end
@@ -111,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]

test/test_tree_2d_eulermulti.jl

@@ -27,13 +27,9 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_2d_dgsem")
       initial_refinement_level = 3,
       tspan = (0.0, 0.001),
       output_directory="out")
-      cmd = string(Base.julia_cmd())
-      coverage = occursin("--code-coverage", cmd) && !occursin("--code-coverage=none", cmd)
-      if !coverage
-        lines = readlines("out/deviations.txt")
-        @test lines[1] == "# iter, simu_time, rho1_min, rho2_min"
-        @test startswith(lines[end], "15")
-      end
+      lines = readlines("out/deviations.txt")
+      @test lines[1] == "# iter, simu_time, rho1_min, rho2_min"
+      @test startswith(lines[end], "1")
   end
 
   @trixi_testset "elixir_eulermulti_ec.jl" begin