Add implementation of subcell bounds check (#1672)

* Add implementation of bounds check

* Shorten code by relocate use of `interval`

* Add note in docstring

* Fix note in docstring of `SubcellLimiterIDP`

* Add test for deviations file

* Adapt last commit to specific setup in test

* Adapt test

* Disable test if `!coverage`

* Remove non-needed code

* Implement first suggestions

* Dispatch `init/finalize_callback()` routines

* Revise docstring

* Rework calculation of bounds

* Fix typo

* Add comment

* Replace construction of `Symbol`s

* Implement suggestions for docstring

examples/tree_2d_dgsem/elixir_euler_shockcapturing_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_positivity.jl

@@ -132,7 +132,9 @@ callbacks           = CallbackSet(summary_callback,
 ###############################################################################
 # run the simulation
 
-stage_callbacks = (SubcellLimiterIDPCorrection(),)
+output_directory = "out"
+stage_callbacks = (SubcellLimiterIDPCorrection(),
+                   BoundsCheckCallback(save_errors=true, interval=100, output_directory=output_directory))
 
 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/Trixi.jl

@@ -232,7 +232,7 @@ export DG,
        SurfaceIntegralUpwind,
        MortarL2
 
-export VolumeIntegralSubcellLimiting,
+export VolumeIntegralSubcellLimiting, BoundsCheckCallback,
        SubcellLimiterIDP, SubcellLimiterIDPCorrection
 
 export nelements, nnodes, nvariables,

src/callbacks_stage/callbacks_stage.jl

@@ -7,6 +7,7 @@
 
 include("positivity_zhang_shu.jl")
 include("subcell_limiter_idp_correction.jl")
+include("subcell_bounds_check.jl")
 # TODO: TrixiShallowWater: move specific limiter file
 include("positivity_shallow_water.jl")
 end # @muladd

src/callbacks_stage/subcell_bounds_check.jl

@@ -0,0 +1,130 @@
+# By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
+# Since these FMAs can increase the performance of many numerical algorithms,
+# we need to opt-in explicitly.
+# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
+@muladd begin
+#! format: noindent
+
+"""
+    BoundsCheckCallback(; output_directory="out", save_errors=false, interval=1)
+
+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.
+For more insights, when setting `save_errors=true` the occurring errors are exported every
+`interval` time steps during the simulation. Then, the maximum deviations since the last
+export are saved in "`output_directory`/deviations.txt".
+The `BoundsCheckCallback` has to be applied as a stage callback for the SSPRK time integration scheme.
+
+!!! note
+    For `SubcellLimiterIDP`, the solution is corrected in the a posteriori correction stage
+    [`SubcellLimiterIDPCorrection`](@ref). So, to check the final solution, this bounds check
+    callback must be called after the correction stage.
+"""
+struct BoundsCheckCallback
+    output_directory::String
+    save_errors::Bool
+    interval::Int
+end
+
+function BoundsCheckCallback(; output_directory = "out", save_errors = false,
+                             interval = 1)
+    BoundsCheckCallback(output_directory, save_errors, interval)
+end
+
+function (callback::BoundsCheckCallback)(u_ode, integrator, stage)
+    mesh, equations, solver, cache = mesh_equations_solver_cache(integrator.p)
+    (; 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)
+    @trixi_timeit timer() "check_bounds" check_bounds(u, mesh, equations, solver, cache,
+                                                      solver.volume_integral, t,
+                                                      iter + 1,
+                                                      callback.output_directory,
+                                                      save_errors)
+end
+
+function check_bounds(u, mesh, equations, solver, cache,
+                      volume_integral::AbstractVolumeIntegral, t, iter,
+                      output_directory, save_errors)
+    return nothing
+end
+
+function check_bounds(u, mesh, equations, solver, cache,
+                      volume_integral::VolumeIntegralSubcellLimiting, t, iter,
+                      output_directory, save_errors)
+    check_bounds(u, mesh, equations, solver, cache, volume_integral.limiter, t, iter,
+                 output_directory, save_errors)
+end
+
+function init_callback(callback::BoundsCheckCallback, semi)
+    init_callback(callback, semi, semi.solver.volume_integral)
+end
+
+init_callback(callback::BoundsCheckCallback, semi, volume_integral::AbstractVolumeIntegral) = nothing
+
+function init_callback(callback::BoundsCheckCallback, semi,
+                       volume_integral::VolumeIntegralSubcellLimiting)
+    init_callback(callback, semi, volume_integral.limiter)
+end
+
+function init_callback(callback::BoundsCheckCallback, semi, limiter::SubcellLimiterIDP)
+    if !callback.save_errors || (callback.interval == 0)
+        return nothing
+    end
+
+    (; 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 positivity
+            for v in limiter.positivity_variables_cons
+                print(f, ", " * string(variables[v]) * "_min")
+            end
+        end
+        println(f)
+    end
+
+    return nothing
+end
+
+function finalize_callback(callback::BoundsCheckCallback, semi)
+    finalize_callback(callback, semi, semi.solver.volume_integral)
+end
+
+finalize_callback(callback::BoundsCheckCallback, semi, volume_integral::AbstractVolumeIntegral) = nothing
+
+function finalize_callback(callback::BoundsCheckCallback, semi,
+                           volume_integral::VolumeIntegralSubcellLimiting)
+    finalize_callback(callback, semi, volume_integral.limiter)
+end
+
+@inline function finalize_callback(callback::BoundsCheckCallback, semi,
+                                   limiter::SubcellLimiterIDP)
+    (; positivity) = limiter
+    (; idp_bounds_delta) = limiter.cache
+    variables = varnames(cons2cons, semi.equations)
+
+    println("─"^100)
+    println("Maximum deviation from bounds:")
+    println("─"^100)
+    if positivity
+        for v in limiter.positivity_variables_cons
+            println(string(variables[v]) * ":\n- positivity: ",
+                    idp_bounds_delta[Symbol(string(v), "_min")][2])
+        end
+    end
+    println("─"^100 * "\n")
+
+    return nothing
+end
+
+include("subcell_bounds_check_2d.jl")
+end # @muladd

src/callbacks_stage/subcell_bounds_check_2d.jl

@@ -0,0 +1,49 @@
+# By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
+# Since these FMAs can increase the performance of many numerical algorithms,
+# we need to opt-in explicitly.
+# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
+@muladd begin
+#! format: noindent
+
+@inline function check_bounds(u, mesh::AbstractMesh{2}, equations, solver, cache,
+                              limiter::SubcellLimiterIDP,
+                              time, iter, output_directory, save_errors)
+    (; positivity) = solver.volume_integral.limiter
+    (; variable_bounds) = limiter.cache.subcell_limiter_coefficients
+    (; idp_bounds_delta) = limiter.cache
+
+    if positivity
+        for v in limiter.positivity_variables_cons
+            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)
+            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
+end
+end # @muladd

src/solvers/dgsem_tree/subcell_limiters.jl

@@ -26,7 +26,7 @@ The bounds are calculated using the low-order FV solution. The positivity limite
 
 !!! note
     This limiter and the correction callback [`SubcellLimiterIDPCorrection`](@ref) only work together.
-    Without the callback, no limiting takes place, leading to a standard flux-differencing DGSEM scheme.
+    Without the callback, no correction takes place, leading to a standard low-order FV scheme.
 
 ## References
 
@@ -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,7 +13,15 @@ function create_cache(limiter::Type{SubcellLimiterIDP}, equations::AbstractEquat
                                                                         nnodes(basis),
                                                                         bound_keys)
 
-    return (; subcell_limiter_coefficients)
+    # 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] = zeros(real(basis), 2)
+    end
+
+    return (; subcell_limiter_coefficients, idp_bounds_delta)
 end
 
 function (limiter::SubcellLimiterIDP)(u::AbstractArray{<:Any, 4}, semi, dg::DGSEM, t,
@@ -60,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

@@ -20,11 +20,16 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_2d_dgsem")
   end
 
   @trixi_testset "elixir_eulermulti_shock_bubble_shockcapturing_subcell_positivity.jl" begin
+    rm("out/deviations.txt", force=true)
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulermulti_shock_bubble_shockcapturing_subcell_positivity.jl"),
       l2   = [81.52845664909304, 2.5455678559421346, 63229.190712645846, 0.19929478404550321, 0.011068604228443425],
       linf = [249.21708417382013, 40.33299887640794, 174205.0118831558, 0.6881458768113586, 0.11274401158173972],
       initial_refinement_level = 3,
-      tspan = (0.0, 0.001))
+      tspan = (0.0, 0.001),
+      output_directory="out")
+      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