Add implementation of subcell bounds check

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

@@ -231,7 +231,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,124 @@
+# 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)
+
+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.
+
+!!! 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, semi)
+    init_callback(callback, semi, semi.solver.volume_integral)
+end
+
+init_callback(callback, semi, volume_integral::AbstractVolumeIntegral) = nothing
+
+function init_callback(callback, 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 index in limiter.positivity_variables_cons
+                print(f, ", $(variables[index])_min")
+            end
+        end
+        println(f)
+    end
+
+    return nothing
+end
+
+function finalize_callback(callback, semi)
+    finalize_callback(callback, semi, semi.solver.volume_integral)
+end
+
+finalize_callback(callback, semi, volume_integral::AbstractVolumeIntegral) = nothing
+
+function finalize_callback(callback, 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("$(variables[v]):\n- positivity: ",
+                    idp_bounds_delta[Symbol("$(v)_min")])
+        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,48 @@
+# 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 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))
+            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
+            end
+        end
+    end
+    if save_errors
+        open("$output_directory/deviations.txt", "a") do f
+            println(f)
+        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
 

src/solvers/dgsem_tree/subcell_limiters_2d.jl

@@ -13,7 +13,12 @@ function create_cache(limiter::Type{SubcellLimiterIDP}, equations::AbstractEquat
                                                                         nnodes(basis),
                                                                         bound_keys)
 
-    return (; subcell_limiter_coefficients)
+    idp_bounds_delta = Dict{Symbol, real(basis)}()
+    for key in bound_keys
+        idp_bounds_delta[key] = zero(real(basis))
+    end
+
+    return (; subcell_limiter_coefficients, idp_bounds_delta)
 end
 
 function (limiter::SubcellLimiterIDP)(u::AbstractArray{<:Any, 4}, semi, dg::DGSEM, t,

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