From 20162dc037bed51b4cea69992eddaae06455251f Mon Sep 17 00:00:00 2001 From: bennibolm Date: Mon, 16 Oct 2023 17:57:29 +0200 Subject: [PATCH 01/17] Add implementation of bounds check --- .../elixir_euler_shockcapturing_subcell.jl | 2 +- ...ubble_shockcapturing_subcell_positivity.jl | 2 +- src/Trixi.jl | 2 +- src/callbacks_stage/callbacks_stage.jl | 1 + src/callbacks_stage/subcell_bounds_check.jl | 123 ++++++++++++++++++ .../subcell_bounds_check_2d.jl | 49 +++++++ src/solvers/dgsem_tree/subcell_limiters_2d.jl | 7 +- 7 files changed, 182 insertions(+), 4 deletions(-) create mode 100644 src/callbacks_stage/subcell_bounds_check.jl create mode 100644 src/callbacks_stage/subcell_bounds_check_2d.jl diff --git a/examples/tree_2d_dgsem/elixir_euler_shockcapturing_subcell.jl b/examples/tree_2d_dgsem/elixir_euler_shockcapturing_subcell.jl index 6b69e4db563..c696e2de416 100644 --- a/examples/tree_2d_dgsem/elixir_euler_shockcapturing_subcell.jl +++ b/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 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 a67eaeb5b2b..a00b75b22e5 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 @@ -132,7 +132,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 diff --git a/src/Trixi.jl b/src/Trixi.jl index b65d03e7975..6414bcc3c42 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -231,7 +231,7 @@ export DG, SurfaceIntegralUpwind, MortarL2 -export VolumeIntegralSubcellLimiting, +export VolumeIntegralSubcellLimiting, BoundsCheckCallback, SubcellLimiterIDP, SubcellLimiterIDPCorrection export nelements, nnodes, nvariables, diff --git a/src/callbacks_stage/callbacks_stage.jl b/src/callbacks_stage/callbacks_stage.jl index 976af327e6f..70d60de7914 100644 --- a/src/callbacks_stage/callbacks_stage.jl +++ b/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 diff --git a/src/callbacks_stage/subcell_bounds_check.jl b/src/callbacks_stage/subcell_bounds_check.jl new file mode 100644 index 00000000000..04b8b4a0652 --- /dev/null +++ b/src/callbacks_stage/subcell_bounds_check.jl @@ -0,0 +1,123 @@ +# 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. +""" +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)) + @trixi_timeit timer() "check_bounds" check_bounds(u, mesh, equations, solver, cache, + t, iter + 1, + callback.output_directory, + save_errors_, callback.interval) +end + +function check_bounds(u, mesh, equations, solver, cache, t, iter, output_directory, + save_errors, interval) + check_bounds(u, mesh, equations, solver, cache, solver.volume_integral, t, iter, + output_directory, save_errors, interval) +end + +function check_bounds(u, mesh, equations, solver, cache, + volume_integral::AbstractVolumeIntegral, + t, iter, output_directory, save_errors, interval) + return nothing +end + +function check_bounds(u, mesh, equations, solver, cache, + volume_integral::VolumeIntegralSubcellLimiting, + t, iter, output_directory, save_errors, interval) + check_bounds(u, mesh, equations, solver, cache, volume_integral.limiter, t, iter, + output_directory, save_errors, interval) +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 diff --git a/src/callbacks_stage/subcell_bounds_check_2d.jl b/src/callbacks_stage/subcell_bounds_check_2d.jl new file mode 100644 index 00000000000..b021f4f2132 --- /dev/null +++ b/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, interval) + (; positivity) = solver.volume_integral.limiter + (; variable_bounds) = limiter.cache.subcell_limiter_coefficients + (; idp_bounds_delta) = limiter.cache + + save_errors_ = save_errors && (iter % interval == 0) + 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 diff --git a/src/solvers/dgsem_tree/subcell_limiters_2d.jl b/src/solvers/dgsem_tree/subcell_limiters_2d.jl index 4497217fb56..690646e9b5a 100644 --- a/src/solvers/dgsem_tree/subcell_limiters_2d.jl +++ b/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, From 3c3c3dab42ffba197aedb0fa5ed9e084adad946f Mon Sep 17 00:00:00 2001 From: bennibolm Date: Tue, 17 Oct 2023 10:08:31 +0200 Subject: [PATCH 02/17] Shorten code by relocate use of `interval` --- src/callbacks_stage/subcell_bounds_check.jl | 18 +++++++++--------- src/callbacks_stage/subcell_bounds_check_2d.jl | 9 ++++----- 2 files changed, 13 insertions(+), 14 deletions(-) diff --git a/src/callbacks_stage/subcell_bounds_check.jl b/src/callbacks_stage/subcell_bounds_check.jl index 04b8b4a0652..919cac3ae23 100644 --- a/src/callbacks_stage/subcell_bounds_check.jl +++ b/src/callbacks_stage/subcell_bounds_check.jl @@ -29,30 +29,30 @@ function (callback::BoundsCheckCallback)(u_ode, integrator, stage) 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) && (stage == length(alg.c)) @trixi_timeit timer() "check_bounds" check_bounds(u, mesh, equations, solver, cache, t, iter + 1, callback.output_directory, - save_errors_, callback.interval) + save_errors_) end function check_bounds(u, mesh, equations, solver, cache, t, iter, output_directory, - save_errors, interval) + save_errors) check_bounds(u, mesh, equations, solver, cache, solver.volume_integral, t, iter, - output_directory, save_errors, interval) + output_directory, save_errors) end function check_bounds(u, mesh, equations, solver, cache, - volume_integral::AbstractVolumeIntegral, - t, iter, output_directory, save_errors, interval) + 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, interval) + 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, interval) + output_directory, save_errors) end function init_callback(callback, semi) diff --git a/src/callbacks_stage/subcell_bounds_check_2d.jl b/src/callbacks_stage/subcell_bounds_check_2d.jl index b021f4f2132..79b0925c444 100644 --- a/src/callbacks_stage/subcell_bounds_check_2d.jl +++ b/src/callbacks_stage/subcell_bounds_check_2d.jl @@ -7,13 +7,12 @@ @inline function check_bounds(u, mesh::AbstractMesh{2}, equations, solver, cache, limiter::SubcellLimiterIDP, - time, iter, output_directory, save_errors, interval) + time, iter, output_directory, save_errors) (; positivity) = solver.volume_integral.limiter (; variable_bounds) = limiter.cache.subcell_limiter_coefficients (; idp_bounds_delta) = limiter.cache - save_errors_ = save_errors && (iter % interval == 0) - if save_errors_ + if save_errors open("$output_directory/deviations.txt", "a") do f print(f, iter, ", ", time) end @@ -30,7 +29,7 @@ variable_bounds[key][i, j, element] - var) end idp_bounds_delta[key] = max(idp_bounds_delta[key], deviation_min) - if save_errors_ + if save_errors deviation_min_ = deviation_min open("$output_directory/deviations.txt", "a") do f print(f, ", ", deviation_min_) @@ -38,7 +37,7 @@ end end end - if save_errors_ + if save_errors open("$output_directory/deviations.txt", "a") do f println(f) end From 21ba6ff95249ce7ca6bd1f64e1d678774f433e31 Mon Sep 17 00:00:00 2001 From: bennibolm Date: Tue, 17 Oct 2023 10:27:35 +0200 Subject: [PATCH 03/17] Add note in docstring --- src/callbacks_stage/subcell_bounds_check.jl | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/src/callbacks_stage/subcell_bounds_check.jl b/src/callbacks_stage/subcell_bounds_check.jl index 919cac3ae23..fb4fa116d02 100644 --- a/src/callbacks_stage/subcell_bounds_check.jl +++ b/src/callbacks_stage/subcell_bounds_check.jl @@ -11,6 +11,11 @@ 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 From 021f96a7c41f08c397ac789a628e60c9b3396563 Mon Sep 17 00:00:00 2001 From: bennibolm Date: Tue, 17 Oct 2023 11:11:48 +0200 Subject: [PATCH 04/17] Fix note in docstring of `SubcellLimiterIDP` --- src/solvers/dgsem_tree/subcell_limiters.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/solvers/dgsem_tree/subcell_limiters.jl b/src/solvers/dgsem_tree/subcell_limiters.jl index 4d9eec25a89..ab519ff8db0 100644 --- a/src/solvers/dgsem_tree/subcell_limiters.jl +++ b/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 From b6d2416bbc296afd32dca2b0997e21a16c1071b7 Mon Sep 17 00:00:00 2001 From: bennibolm Date: Tue, 17 Oct 2023 15:09:39 +0200 Subject: [PATCH 05/17] Add test for deviations file --- ...rmulti_shock_bubble_shockcapturing_subcell_positivity.jl | 4 +++- src/callbacks_stage/subcell_bounds_check.jl | 3 ++- test/test_tree_2d_eulermulti.jl | 6 +++++- 3 files changed, 10 insertions(+), 3 deletions(-) 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 a00b75b22e5..c5a7a5932e6 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 @@ -132,7 +132,9 @@ callbacks = CallbackSet(summary_callback, ############################################################################### # run the simulation -stage_callbacks = (SubcellLimiterIDPCorrection(), BoundsCheckCallback(save_errors=false)) +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 diff --git a/src/callbacks_stage/subcell_bounds_check.jl b/src/callbacks_stage/subcell_bounds_check.jl index fb4fa116d02..791b6d2a4c3 100644 --- a/src/callbacks_stage/subcell_bounds_check.jl +++ b/src/callbacks_stage/subcell_bounds_check.jl @@ -34,7 +34,8 @@ function (callback::BoundsCheckCallback)(u_ode, integrator, stage) u = wrap_array(u_ode, mesh, equations, solver, cache) save_errors_ = callback.save_errors && (callback.interval > 0) && - (iter % callback.interval == 0) && (stage == length(alg.c)) + (stage == length(alg.c)) && + (iter % callback.interval == 0 || integrator.finalstep) @trixi_timeit timer() "check_bounds" check_bounds(u, mesh, equations, solver, cache, t, iter + 1, callback.output_directory, diff --git a/test/test_tree_2d_eulermulti.jl b/test/test_tree_2d_eulermulti.jl index 606afca1034..eafd9205b98 100644 --- a/test/test_tree_2d_eulermulti.jl +++ b/test/test_tree_2d_eulermulti.jl @@ -24,7 +24,11 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_2d_dgsem") 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") + file = "out/deviations.txt" + lines = readlines(file) + @assert lines[end] == "139, 0.00991608811086362, 0.0, 3.122502256758253e-17" end @trixi_testset "elixir_eulermulti_ec.jl" begin From 5249783b9c136fa10eeabb35378b669961eef130 Mon Sep 17 00:00:00 2001 From: bennibolm Date: Tue, 17 Oct 2023 15:11:36 +0200 Subject: [PATCH 06/17] Adapt last commit to specific setup in test --- test/test_tree_2d_eulermulti.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_tree_2d_eulermulti.jl b/test/test_tree_2d_eulermulti.jl index eafd9205b98..7be83b8315e 100644 --- a/test/test_tree_2d_eulermulti.jl +++ b/test/test_tree_2d_eulermulti.jl @@ -28,7 +28,7 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_2d_dgsem") output_directory="out") file = "out/deviations.txt" lines = readlines(file) - @assert lines[end] == "139, 0.00991608811086362, 0.0, 3.122502256758253e-17" + @assert lines[end] == "15, 0.0009595796113231045, 0.0, 0.0" end @trixi_testset "elixir_eulermulti_ec.jl" begin From ebf37af9d84007fa605c9453bebb8fb154b9e8f1 Mon Sep 17 00:00:00 2001 From: bennibolm Date: Tue, 17 Oct 2023 15:41:03 +0200 Subject: [PATCH 07/17] Adapt test --- test/test_tree_2d_eulermulti.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/test/test_tree_2d_eulermulti.jl b/test/test_tree_2d_eulermulti.jl index 7be83b8315e..9cb90121723 100644 --- a/test/test_tree_2d_eulermulti.jl +++ b/test/test_tree_2d_eulermulti.jl @@ -20,6 +20,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) @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], @@ -28,7 +29,7 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_2d_dgsem") output_directory="out") file = "out/deviations.txt" lines = readlines(file) - @assert lines[end] == "15, 0.0009595796113231045, 0.0, 0.0" + @assert startswith(lines[end], "15") end @trixi_testset "elixir_eulermulti_ec.jl" begin From 6c7d637747c14e53039370943c39723a68fe043f Mon Sep 17 00:00:00 2001 From: bennibolm Date: Wed, 18 Oct 2023 10:21:06 +0200 Subject: [PATCH 08/17] Disable test if `!coverage` --- test/test_tree_2d_eulermulti.jl | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/test/test_tree_2d_eulermulti.jl b/test/test_tree_2d_eulermulti.jl index 9cb90121723..126532aa0ee 100644 --- a/test/test_tree_2d_eulermulti.jl +++ b/test/test_tree_2d_eulermulti.jl @@ -27,9 +27,13 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_2d_dgsem") initial_refinement_level = 3, tspan = (0.0, 0.001), output_directory="out") - file = "out/deviations.txt" - lines = readlines(file) - @assert startswith(lines[end], "15") + 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 end @trixi_testset "elixir_eulermulti_ec.jl" begin From 97db537f8ed19961765bd3d20fc367ac3251cc62 Mon Sep 17 00:00:00 2001 From: bennibolm Date: Wed, 18 Oct 2023 12:15:41 +0200 Subject: [PATCH 09/17] Remove non-needed code --- src/callbacks_stage/subcell_bounds_check.jl | 9 ++------- 1 file changed, 2 insertions(+), 7 deletions(-) diff --git a/src/callbacks_stage/subcell_bounds_check.jl b/src/callbacks_stage/subcell_bounds_check.jl index 791b6d2a4c3..581a9cbb82b 100644 --- a/src/callbacks_stage/subcell_bounds_check.jl +++ b/src/callbacks_stage/subcell_bounds_check.jl @@ -37,17 +37,12 @@ function (callback::BoundsCheckCallback)(u_ode, integrator, stage) (stage == length(alg.c)) && (iter % callback.interval == 0 || integrator.finalstep) @trixi_timeit timer() "check_bounds" check_bounds(u, mesh, equations, solver, cache, - t, iter + 1, + solver.volume_integral, t, + iter + 1, callback.output_directory, save_errors_) end -function check_bounds(u, mesh, equations, solver, cache, t, iter, output_directory, - save_errors) - check_bounds(u, mesh, equations, solver, cache, solver.volume_integral, t, iter, - output_directory, save_errors) -end - function check_bounds(u, mesh, equations, solver, cache, volume_integral::AbstractVolumeIntegral, t, iter, output_directory, save_errors) From 4be4d61c798a2521ceda5c44b94882180358a109 Mon Sep 17 00:00:00 2001 From: bennibolm Date: Wed, 18 Oct 2023 15:43:08 +0200 Subject: [PATCH 10/17] Implement first suggestions --- src/callbacks_stage/subcell_bounds_check.jl | 8 ++++---- test/test_tree_2d_eulermulti.jl | 10 +++------- 2 files changed, 7 insertions(+), 11 deletions(-) diff --git a/src/callbacks_stage/subcell_bounds_check.jl b/src/callbacks_stage/subcell_bounds_check.jl index 581a9cbb82b..1b605d19428 100644 --- a/src/callbacks_stage/subcell_bounds_check.jl +++ b/src/callbacks_stage/subcell_bounds_check.jl @@ -33,14 +33,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, diff --git a/test/test_tree_2d_eulermulti.jl b/test/test_tree_2d_eulermulti.jl index 126532aa0ee..0ae46a92ef8 100644 --- a/test/test_tree_2d_eulermulti.jl +++ b/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 From c9b83a6eb4f9308d3261ac3b8055c3c08e851402 Mon Sep 17 00:00:00 2001 From: bennibolm Date: Thu, 19 Oct 2023 11:39:13 +0200 Subject: [PATCH 11/17] Dispatch `init/finalize_callback()` routines --- src/callbacks_stage/subcell_bounds_check.jl | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/src/callbacks_stage/subcell_bounds_check.jl b/src/callbacks_stage/subcell_bounds_check.jl index 1b605d19428..8be1c8b3800 100644 --- a/src/callbacks_stage/subcell_bounds_check.jl +++ b/src/callbacks_stage/subcell_bounds_check.jl @@ -56,13 +56,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 @@ -89,13 +90,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 From 5f01e42eb8498f4262983004348289bae0ed775d Mon Sep 17 00:00:00 2001 From: bennibolm Date: Thu, 19 Oct 2023 14:58:25 +0200 Subject: [PATCH 12/17] Revise docstring --- src/callbacks_stage/subcell_bounds_check.jl | 15 ++++++++++----- 1 file changed, 10 insertions(+), 5 deletions(-) diff --git a/src/callbacks_stage/subcell_bounds_check.jl b/src/callbacks_stage/subcell_bounds_check.jl index 8be1c8b3800..81e7007f5a3 100644 --- a/src/callbacks_stage/subcell_bounds_check.jl +++ b/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 occuring 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 @@ -80,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, ", $(variables[v])_min") end end println(f) From 9e9e99e1c591f462f184795503d74796211bdfe0 Mon Sep 17 00:00:00 2001 From: bennibolm Date: Thu, 19 Oct 2023 14:59:39 +0200 Subject: [PATCH 13/17] Rework calculation of bounds --- src/callbacks_stage/subcell_bounds_check.jl | 2 +- .../subcell_bounds_check_2d.jl | 31 ++++++++++--------- src/solvers/dgsem_tree/subcell_limiters_2d.jl | 4 +-- 3 files changed, 19 insertions(+), 18 deletions(-) diff --git a/src/callbacks_stage/subcell_bounds_check.jl b/src/callbacks_stage/subcell_bounds_check.jl index 81e7007f5a3..d62cf501054 100644 --- a/src/callbacks_stage/subcell_bounds_check.jl +++ b/src/callbacks_stage/subcell_bounds_check.jl @@ -118,7 +118,7 @@ end if positivity for v in limiter.positivity_variables_cons println("$(variables[v]):\n- positivity: ", - idp_bounds_delta[Symbol("$(v)_min")]) + idp_bounds_delta[Symbol("$(v)_min")][2]) end end println("─"^100 * "\n") diff --git a/src/callbacks_stage/subcell_bounds_check_2d.jl b/src/callbacks_stage/subcell_bounds_check_2d.jl index 79b0925c444..2df10d0a757 100644 --- a/src/callbacks_stage/subcell_bounds_check_2d.jl +++ b/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)) + 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("$(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 diff --git a/src/solvers/dgsem_tree/subcell_limiters_2d.jl b/src/solvers/dgsem_tree/subcell_limiters_2d.jl index 690646e9b5a..4d32272898b 100644 --- a/src/solvers/dgsem_tree/subcell_limiters_2d.jl +++ b/src/solvers/dgsem_tree/subcell_limiters_2d.jl @@ -13,9 +13,9 @@ function create_cache(limiter::Type{SubcellLimiterIDP}, equations::AbstractEquat nnodes(basis), bound_keys) - idp_bounds_delta = Dict{Symbol, real(basis)}() + 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) From a27ac02e7467387338293f3264ec0fb7d3fe0256 Mon Sep 17 00:00:00 2001 From: bennibolm Date: Thu, 19 Oct 2023 15:02:06 +0200 Subject: [PATCH 14/17] Fix typo --- src/callbacks_stage/subcell_bounds_check.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/callbacks_stage/subcell_bounds_check.jl b/src/callbacks_stage/subcell_bounds_check.jl index 81e7007f5a3..66e11dc9222 100644 --- a/src/callbacks_stage/subcell_bounds_check.jl +++ b/src/callbacks_stage/subcell_bounds_check.jl @@ -12,7 +12,7 @@ Subcell limiting techniques with [`SubcellLimiterIDP`](@ref) are constructed to 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 occuring errors can be exported with `save_errors==true` +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. From aa44953a52df0be9fbf93017035778943477e8a2 Mon Sep 17 00:00:00 2001 From: bennibolm Date: Thu, 19 Oct 2023 15:26:01 +0200 Subject: [PATCH 15/17] Add comment --- src/solvers/dgsem_tree/subcell_limiters_2d.jl | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/solvers/dgsem_tree/subcell_limiters_2d.jl b/src/solvers/dgsem_tree/subcell_limiters_2d.jl index 4d32272898b..e529394c660 100644 --- a/src/solvers/dgsem_tree/subcell_limiters_2d.jl +++ b/src/solvers/dgsem_tree/subcell_limiters_2d.jl @@ -13,6 +13,9 @@ function create_cache(limiter::Type{SubcellLimiterIDP}, equations::AbstractEquat nnodes(basis), bound_keys) + # 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) From c509ef5a24806aba20c78af91ae57cb9d63ee766 Mon Sep 17 00:00:00 2001 From: bennibolm Date: Fri, 20 Oct 2023 12:50:22 +0200 Subject: [PATCH 16/17] Replace construction of `Symbol`s --- src/callbacks_stage/subcell_bounds_check.jl | 6 +++--- src/callbacks_stage/subcell_bounds_check_2d.jl | 4 ++-- src/solvers/dgsem_tree/subcell_limiters.jl | 2 +- src/solvers/dgsem_tree/subcell_limiters_2d.jl | 2 +- 4 files changed, 7 insertions(+), 7 deletions(-) diff --git a/src/callbacks_stage/subcell_bounds_check.jl b/src/callbacks_stage/subcell_bounds_check.jl index d076a4aaa0f..55b5841c227 100644 --- a/src/callbacks_stage/subcell_bounds_check.jl +++ b/src/callbacks_stage/subcell_bounds_check.jl @@ -86,7 +86,7 @@ function init_callback(callback::BoundsCheckCallback, semi, limiter::SubcellLimi print(f, "# iter, simu_time") if positivity for v in limiter.positivity_variables_cons - print(f, ", $(variables[v])_min") + print(f, ", " * string(variables[v]) * "_min") end end println(f) @@ -117,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")][2]) + println(string(variables[v]) * ":\n- positivity: ", + idp_bounds_delta[Symbol(string(v), "_min")][2]) end end println("─"^100 * "\n") diff --git a/src/callbacks_stage/subcell_bounds_check_2d.jl b/src/callbacks_stage/subcell_bounds_check_2d.jl index 2df10d0a757..8159becb503 100644 --- a/src/callbacks_stage/subcell_bounds_check_2d.jl +++ b/src/callbacks_stage/subcell_bounds_check_2d.jl @@ -14,7 +14,7 @@ if positivity for v in limiter.positivity_variables_cons - key = Symbol("$(v)_min") + key = Symbol(string(v), "_min") deviation = idp_bounds_delta[key] for element in eachelement(solver, cache), j in eachnode(solver), i in eachnode(solver) @@ -32,7 +32,7 @@ print(f, iter, ", ", time) if positivity for v in limiter.positivity_variables_cons - key = Symbol("$(v)_min") + key = Symbol(string(v), "_min") print(f, ", ", idp_bounds_delta[key][1]) end end diff --git a/src/solvers/dgsem_tree/subcell_limiters.jl b/src/solvers/dgsem_tree/subcell_limiters.jl index ab519ff8db0..3a508dd180e 100644 --- a/src/solvers/dgsem_tree/subcell_limiters.jl +++ b/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) diff --git a/src/solvers/dgsem_tree/subcell_limiters_2d.jl b/src/solvers/dgsem_tree/subcell_limiters_2d.jl index e529394c660..5e00ab4e903 100644 --- a/src/solvers/dgsem_tree/subcell_limiters_2d.jl +++ b/src/solvers/dgsem_tree/subcell_limiters_2d.jl @@ -68,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] From f4f4146fd644068041e376acd3b482bad656b26e Mon Sep 17 00:00:00 2001 From: bennibolm Date: Fri, 20 Oct 2023 14:43:53 +0200 Subject: [PATCH 17/17] Implement suggestions for docstring --- src/callbacks_stage/subcell_bounds_check.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/callbacks_stage/subcell_bounds_check.jl b/src/callbacks_stage/subcell_bounds_check.jl index 55b5841c227..c86f266147c 100644 --- a/src/callbacks_stage/subcell_bounds_check.jl +++ b/src/callbacks_stage/subcell_bounds_check.jl @@ -12,10 +12,10 @@ Subcell limiting techniques with [`SubcellLimiterIDP`](@ref) are constructed to 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 +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". -It has to be implied as a stage callback for SSPRK. +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