Outsource saving of bounds check errors (#1824)

* Outsource saving of bounds check errors

* Increase coverage by activate saving of bounds check errors for local limiting

* Fix test for saving errors with coverage

* Remove empty line

* Implement suggestions

* Use `joinpath`

* Add test with coverage

* Fix I/O test for modified elixir

* Add more accurate test

* Redo last change; Default maxiters with coverage is 1

* Change order of lines

* Save deviations if maxiters is reached

* Add comment

* Save deviations at correct time steps

* Add test for saving deviations for positivity limiter (non-linear)

* Add comment

---------

Co-authored-by: Michael Schlottke-Lakemper <michael@sloede.com>

docs/literate/src/files/subcell_shock_capturing.jl

@@ -277,11 +277,10 @@ plot(sol)
 # timestep and simulation time.
 # ````
 # iter, simu_time, rho_min, rho_max
-# 1, 0.0, 0.0, 0.0
-# 101, 0.29394033217556337, 0.0, 0.0
-# 201, 0.6012597465597065, 0.0, 0.0
-# 301, 0.9559096690030839, 0.0, 0.0
-# 401, 1.3674274981949077, 0.0, 0.0
-# 501, 1.8395301696603052, 0.0, 0.0
+# 100, 0.29103427131404924, 0.0, 0.0
+# 200, 0.5980281923063808, 0.0, 0.0
+# 300, 0.9520853560765293, 0.0, 0.0
+# 400, 1.3630295622683186, 0.0, 0.0
+# 500, 1.8344999624013498, 0.0, 0.0
 # 532, 1.9974179806990118, 0.0, 0.0
 # ````

examples/tree_2d_dgsem/elixir_euler_kelvin_helmholtz_instability_sc_subcell.jl

@@ -83,7 +83,9 @@ callbacks = CallbackSet(summary_callback,
 ###############################################################################
 # run the simulation
 
-stage_callbacks = (SubcellLimiterIDPCorrection(), BoundsCheckCallback(save_errors = false))
+stage_callbacks = (SubcellLimiterIDPCorrection(),
+                   BoundsCheckCallback(save_errors = false, interval = 100))
+# `interval` is used when calling this elixir in the tests with `save_errors=true`.
 
 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

@@ -85,7 +85,9 @@ callbacks = CallbackSet(summary_callback,
 ###############################################################################
 # run the simulation
 
-stage_callbacks = (SubcellLimiterIDPCorrection(), BoundsCheckCallback(save_errors = false))
+stage_callbacks = (SubcellLimiterIDPCorrection(),
+                   BoundsCheckCallback(save_errors = false, interval = 100))
+# `interval` is used when calling this elixir in the tests with `save_errors=true`.
 
 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

@@ -133,10 +133,9 @@ callbacks = CallbackSet(summary_callback,
 ###############################################################################
 # run the simulation
 
-output_directory = "out"
 stage_callbacks = (SubcellLimiterIDPCorrection(),
-                   BoundsCheckCallback(save_errors = true, interval = 100,
-                                       output_directory = output_directory))
+                   BoundsCheckCallback(save_errors = false, interval = 100))
+# `interval` is used when calling this elixir in the tests with `save_errors=true`.
 
 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

@@ -38,35 +38,37 @@ function (callback::BoundsCheckCallback)(u_ode, integrator, stage)
     (; t, iter, alg) = integrator
     u = wrap_array(u_ode, mesh, equations, solver, cache)
 
+    @trixi_timeit timer() "check_bounds" check_bounds(u, mesh, equations, solver, cache,
+                                                      solver.volume_integral)
+
     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)
+                  ((iter + 1) % callback.interval == 0 ||   # Every `interval` time steps
+                   integrator.finalstep ||                  # Planned last time step
+                   (iter + 1) >= integrator.opts.maxiters)  # Maximum iterations reached
+    if save_errors
+        @trixi_timeit timer() "save_errors" save_bounds_check_errors(callback.output_directory,
+                                                                     u, t, iter + 1,
+                                                                     equations,
+                                                                     solver.volume_integral)
+    end
 end
 
-function check_bounds(u, mesh, equations, solver, cache,
-                      volume_integral::AbstractVolumeIntegral, t, iter,
-                      output_directory, save_errors)
-    return nothing
+@inline function check_bounds(u, mesh, equations, solver, cache,
+                              volume_integral::VolumeIntegralSubcellLimiting)
+    check_bounds(u, mesh, equations, solver, cache, volume_integral.limiter)
 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)
+@inline function save_bounds_check_errors(output_directory, u, t, iter, equations,
+                                          volume_integral::VolumeIntegralSubcellLimiting)
+    save_bounds_check_errors(output_directory, u, t, iter, equations,
+                             volume_integral.limiter)
 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)
@@ -116,8 +118,6 @@ 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)

src/callbacks_stage/subcell_bounds_check_2d.jl

@@ -6,8 +6,7 @@
 #! format: noindent
 
 @inline function check_bounds(u, mesh::AbstractMesh{2}, equations, solver, cache,
-                              limiter::SubcellLimiterIDP,
-                              time, iter, output_directory, save_errors)
+                              limiter::SubcellLimiterIDP)
     (; local_twosided, positivity, local_onesided) = solver.volume_integral.limiter
     (; variable_bounds) = limiter.cache.subcell_limiter_coefficients
     (; idp_bounds_delta_local, idp_bounds_delta_global) = limiter.cache
@@ -101,42 +100,47 @@
                                            idp_bounds_delta_local[key])
     end
 
-    if save_errors
-        # Print to output file
-        open("$output_directory/deviations.txt", "a") do f
-            print(f, iter, ", ", time)
-            if local_twosided
-                for v in limiter.local_twosided_variables_cons
-                    v_string = string(v)
-                    print(f, ", ", idp_bounds_delta_local[Symbol(v_string, "_min")],
-                          ", ", idp_bounds_delta_local[Symbol(v_string, "_max")])
-                end
+    return nothing
+end
+
+@inline function save_bounds_check_errors(output_directory, u, time, iter, equations,
+                                          limiter::SubcellLimiterIDP)
+    (; local_twosided, positivity, local_onesided) = limiter
+    (; idp_bounds_delta_local) = limiter.cache
+
+    # Print to output file
+    open(joinpath(output_directory, "deviations.txt"), "a") do f
+        print(f, iter, ", ", time)
+        if local_twosided
+            for v in limiter.local_twosided_variables_cons
+                v_string = string(v)
+                print(f, ", ", idp_bounds_delta_local[Symbol(v_string, "_min")],
+                      ", ", idp_bounds_delta_local[Symbol(v_string, "_max")])
             end
-            if local_onesided
-                for (variable, min_or_max) in limiter.local_onesided_variables_nonlinear
-                    print(f, ", ",
-                          idp_bounds_delta_local[Symbol(string(variable), "_",
-                                                        string(min_or_max))][stride_size])
-                end
+        end
+        if local_onesided
+            for (variable, min_or_max) in limiter.local_onesided_variables_nonlinear
+                key = Symbol(string(variable), "_", string(min_or_max))
+                print(f, ", ", idp_bounds_delta_local[key])
             end
-            if positivity
-                for v in limiter.positivity_variables_cons
-                    if v in limiter.local_twosided_variables_cons
-                        continue
-                    end
-                    print(f, ", ", idp_bounds_delta_local[Symbol(string(v), "_min")])
-                end
-                for variable in limiter.positivity_variables_nonlinear
-                    print(f, ", ",
-                          idp_bounds_delta_local[Symbol(string(variable), "_min")])
+        end
+        if positivity
+            for v in limiter.positivity_variables_cons
+                if v in limiter.local_twosided_variables_cons
+                    continue
                 end
+                print(f, ", ", idp_bounds_delta_local[Symbol(string(v), "_min")])
+            end
+            for variable in limiter.positivity_variables_nonlinear
+                print(f, ", ", idp_bounds_delta_local[Symbol(string(variable), "_min")])
             end
-            println(f)
-        end
-        # Reset local maximum deviations
-        for (key, _) in idp_bounds_delta_local
-            idp_bounds_delta_local[key] = zero(eltype(idp_bounds_delta_local[key]))
         end
+        println(f)
+    end
+
+    # Reset local maximum deviations
+    for (key, _) in idp_bounds_delta_local
+        idp_bounds_delta_local[key] = zero(eltype(idp_bounds_delta_local[key]))
     end
 
     return nothing

test/test_tree_2d_euler.jl

@@ -400,6 +400,7 @@ end
 end
 
 @trixi_testset "elixir_euler_sedov_blast_wave_sc_subcell.jl" begin
+    rm(joinpath("out", "deviations.txt"), force = true)
     @test_trixi_include(joinpath(EXAMPLES_DIR,
                                  "elixir_euler_sedov_blast_wave_sc_subcell.jl"),
                         l2=[
@@ -416,7 +417,20 @@ end
                         ],
                         tspan=(0.0, 1.0),
                         initial_refinement_level=4,
-                        coverage_override=(maxiters = 6,))
+                        coverage_override=(maxiters = 6,),
+                        save_errors=true)
+    lines = readlines(joinpath("out", "deviations.txt"))
+    @test lines[1] == "# iter, simu_time, rho_min, rho_max, entropy_guermond_etal_min"
+    cmd = string(Base.julia_cmd())
+    coverage = occursin("--code-coverage", cmd) &&
+               !occursin("--code-coverage=none", cmd)
+    if coverage
+        # Run with coverage takes 6 time steps.
+        @test startswith(lines[end], "6")
+    else
+        # Run without coverage takes 89 time steps.
+        @test startswith(lines[end], "89")
+    end
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)
     let
@@ -614,6 +628,7 @@ end
 end
 
 @trixi_testset "elixir_euler_kelvin_helmholtz_instability_sc_subcell.jl" begin
+    rm(joinpath("out", "deviations.txt"), force = true)
     @test_trixi_include(joinpath(EXAMPLES_DIR,
                                  "elixir_euler_kelvin_helmholtz_instability_sc_subcell.jl"),
                         l2=[
@@ -628,7 +643,13 @@ end
                             0.5822547982757897,
                             0.7300051017382696,
                         ],
-                        tspan=(0.0, 2.0))
+                        tspan=(0.0, 2.0),
+                        coverage_override=(maxiters = 7,),
+                        save_errors=true)
+    lines = readlines(joinpath("out", "deviations.txt"))
+    @test lines[1] == "# iter, simu_time, rho_min, pressure_min"
+    # Run without (with) coverage takes 745 (7) time steps
+    @test startswith(lines[end], "7")
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)
     let

test/test_tree_2d_eulermulti.jl

@@ -61,7 +61,7 @@ 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)
+        rm(joinpath("out", "deviations.txt"), force = true)
         @test_trixi_include(joinpath(EXAMPLES_DIR,
                                      "elixir_eulermulti_shock_bubble_shockcapturing_subcell_positivity.jl"),
                             l2=[
@@ -80,9 +80,10 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_2d_dgsem")
                             ],
                             initial_refinement_level=3,
                             tspan=(0.0, 0.001),
-                            output_directory="out")
-        lines = readlines("out/deviations.txt")
+                            save_errors=true)
+        lines = readlines(joinpath("out", "deviations.txt"))
         @test lines[1] == "# iter, simu_time, rho1_min, rho2_min"
+        # Runs with and without coverage take 1 and 15 time steps.
         @test startswith(lines[end], "1")
         # Ensure that we do not have excessive memory allocations
         # (e.g., from type instabilities)