diff --git a/docs/literate/src/files/subcell_shock_capturing.jl b/docs/literate/src/files/subcell_shock_capturing.jl index 8a98fdae28..b4d0896536 100644 --- a/docs/literate/src/files/subcell_shock_capturing.jl +++ b/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 # ```` diff --git a/examples/tree_2d_dgsem/elixir_euler_kelvin_helmholtz_instability_sc_subcell.jl b/examples/tree_2d_dgsem/elixir_euler_kelvin_helmholtz_instability_sc_subcell.jl index 1817672778..9e9fb45e7d 100644 --- a/examples/tree_2d_dgsem/elixir_euler_kelvin_helmholtz_instability_sc_subcell.jl +++ b/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 diff --git a/examples/tree_2d_dgsem/elixir_euler_sedov_blast_wave_sc_subcell.jl b/examples/tree_2d_dgsem/elixir_euler_sedov_blast_wave_sc_subcell.jl index 6cbbe4eb4e..2089d35397 100644 --- a/examples/tree_2d_dgsem/elixir_euler_sedov_blast_wave_sc_subcell.jl +++ b/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 diff --git a/examples/tree_2d_dgsem/elixir_eulermulti_shock_bubble_shockcapturing_subcell_positivity.jl b/examples/tree_2d_dgsem/elixir_eulermulti_shock_bubble_shockcapturing_subcell_positivity.jl index 78ff47e255..be14c448e4 100644 --- a/examples/tree_2d_dgsem/elixir_eulermulti_shock_bubble_shockcapturing_subcell_positivity.jl +++ b/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 diff --git a/src/callbacks_stage/subcell_bounds_check.jl b/src/callbacks_stage/subcell_bounds_check.jl index ba193ab299..3f3e151436 100644 --- a/src/callbacks_stage/subcell_bounds_check.jl +++ b/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) diff --git a/src/callbacks_stage/subcell_bounds_check_2d.jl b/src/callbacks_stage/subcell_bounds_check_2d.jl index 0f713a296e..c81ebc970a 100644 --- a/src/callbacks_stage/subcell_bounds_check_2d.jl +++ b/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 diff --git a/test/test_tree_2d_euler.jl b/test/test_tree_2d_euler.jl index efd4058031..a004d1452b 100644 --- a/test/test_tree_2d_euler.jl +++ b/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 diff --git a/test/test_tree_2d_eulermulti.jl b/test/test_tree_2d_eulermulti.jl index 0aaa9be1c5..5b98461168 100644 --- a/test/test_tree_2d_eulermulti.jl +++ b/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)