From 535f52004f8799c0ae936d33dcdc630a969ced68 Mon Sep 17 00:00:00 2001 From: bennibolm Date: Sun, 22 Oct 2023 19:22:26 +0200 Subject: [PATCH] Add bounds check for local limiting --- .../elixir_euler_blast_wave_sc_subcell.jl | 2 +- ...lixir_euler_sedov_blast_wave_sc_subcell.jl | 2 +- ...ck_bubble_shockcapturing_subcell_minmax.jl | 2 +- src/callbacks_stage/subcell_bounds_check.jl | 24 +++++++++++- .../subcell_bounds_check_2d.jl | 38 +++++++++++++++++-- src/solvers/dgsem_tree/subcell_limiters.jl | 5 ++- 6 files changed, 63 insertions(+), 10 deletions(-) diff --git a/examples/tree_2d_dgsem/elixir_euler_blast_wave_sc_subcell.jl b/examples/tree_2d_dgsem/elixir_euler_blast_wave_sc_subcell.jl index e7b701611ad..792a813dc36 100644 --- a/examples/tree_2d_dgsem/elixir_euler_blast_wave_sc_subcell.jl +++ b/examples/tree_2d_dgsem/elixir_euler_blast_wave_sc_subcell.jl @@ -88,7 +88,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_euler_sedov_blast_wave_sc_subcell.jl b/examples/tree_2d_dgsem/elixir_euler_sedov_blast_wave_sc_subcell.jl index 21de64aee4c..eae24945c65 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 @@ -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_minmax.jl b/examples/tree_2d_dgsem/elixir_eulermulti_shock_bubble_shockcapturing_subcell_minmax.jl index 4184e04de3d..4e5ec4eb8d1 100644 --- a/examples/tree_2d_dgsem/elixir_eulermulti_shock_bubble_shockcapturing_subcell_minmax.jl +++ b/examples/tree_2d_dgsem/elixir_eulermulti_shock_bubble_shockcapturing_subcell_minmax.jl @@ -131,7 +131,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/callbacks_stage/subcell_bounds_check.jl b/src/callbacks_stage/subcell_bounds_check.jl index c86f266147c..d7e30ab1621 100644 --- a/src/callbacks_stage/subcell_bounds_check.jl +++ b/src/callbacks_stage/subcell_bounds_check.jl @@ -77,15 +77,24 @@ function init_callback(callback::BoundsCheckCallback, semi, limiter::SubcellLimi return nothing end - (; positivity) = limiter + (; local_minmax, 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 local_minmax + for v in limiter.local_minmax_variables_cons + variable_string = string(variables[v]) + print(f, ", " * variable_string * "_min, " * variable_string * "_max") + end + end if positivity for v in limiter.positivity_variables_cons + if v in limiter.local_minmax_variables_cons + continue + end print(f, ", " * string(variables[v]) * "_min") end end @@ -108,15 +117,26 @@ end @inline function finalize_callback(callback::BoundsCheckCallback, semi, limiter::SubcellLimiterIDP) - (; positivity) = limiter + (; local_minmax, positivity) = limiter (; idp_bounds_delta) = limiter.cache variables = varnames(cons2cons, semi.equations) println("─"^100) println("Maximum deviation from bounds:") println("─"^100) + if local_minmax + for v in limiter.local_minmax_variables_cons + v_string = string(v) + println("$(variables[v]):") + println("-lower bound: ", idp_bounds_delta[Symbol(v_string, "_min")][2]) + println("-upper bound: ", idp_bounds_delta[Symbol(v_string, "_max")][2]) + end + end if positivity for v in limiter.positivity_variables_cons + if v in limiter.local_minmax_variables_cons + continue + end println(string(variables[v]) * ":\n- positivity: ", idp_bounds_delta[Symbol(string(v), "_min")][2]) end diff --git a/src/callbacks_stage/subcell_bounds_check_2d.jl b/src/callbacks_stage/subcell_bounds_check_2d.jl index 8159becb503..654b4b54499 100644 --- a/src/callbacks_stage/subcell_bounds_check_2d.jl +++ b/src/callbacks_stage/subcell_bounds_check_2d.jl @@ -8,12 +8,35 @@ @inline function check_bounds(u, mesh::AbstractMesh{2}, equations, solver, cache, limiter::SubcellLimiterIDP, time, iter, output_directory, save_errors) - (; positivity) = solver.volume_integral.limiter + (; local_minmax, positivity) = solver.volume_integral.limiter (; variable_bounds) = limiter.cache.subcell_limiter_coefficients (; idp_bounds_delta) = limiter.cache + if local_minmax + for v in limiter.local_minmax_variables_cons + v_string = string(v) + key_min = Symbol(v_string, "_min") + key_max = Symbol(v_string, "_max") + deviation_min = idp_bounds_delta[key_min] + deviation_max = idp_bounds_delta[key_max] + for element in eachelement(solver, cache), j in eachnode(solver), + i in eachnode(solver) + + var = u[v, i, j, element] + deviation_min[1] = max(deviation_min[1], + variable_bounds[key_min][i, j, element] - var) + deviation_max[1] = max(deviation_max[1], + var - variable_bounds[key_max][i, j, element]) + end + deviation_min[2] = max(deviation_min[2], deviation_min[1]) + deviation_max[2] = max(deviation_max[2], deviation_max[1]) + end + end if positivity for v in limiter.positivity_variables_cons + if v in limiter.local_minmax_variables_cons + continue + end key = Symbol(string(v), "_min") deviation = idp_bounds_delta[key] for element in eachelement(solver, cache), j in eachnode(solver), @@ -30,10 +53,19 @@ # Print to output file open("$output_directory/deviations.txt", "a") do f print(f, iter, ", ", time) + if local_minmax + for v in limiter.local_minmax_variables_cons + v_string = string(v) + print(f, ", ", idp_bounds_delta[Symbol(v_string, "_min")][1], + idp_bounds_delta[Symbol(v_string, "_max")][1]) + end + end if positivity for v in limiter.positivity_variables_cons - key = Symbol(string(v), "_min") - print(f, ", ", idp_bounds_delta[key][1]) + if v in limiter.local_minmax_variables_cons + continue + end + print(f, ", ", idp_bounds_delta[Symbol(string(v), "_min")][1]) end end println(f) diff --git a/src/solvers/dgsem_tree/subcell_limiters.jl b/src/solvers/dgsem_tree/subcell_limiters.jl index 2cec0ebf1ef..1c711c9c017 100644 --- a/src/solvers/dgsem_tree/subcell_limiters.jl +++ b/src/solvers/dgsem_tree/subcell_limiters.jl @@ -62,8 +62,9 @@ function SubcellLimiterIDP(equations::AbstractEquations, basis; bound_keys = () if local_minmax for v in local_minmax_variables_cons - bound_keys = (bound_keys..., Symbol(string(v), "_min"), - Symbol(string(v), "_max")) + v_string = string(v) + bound_keys = (bound_keys..., Symbol(v_string, "_min"), + Symbol(v_string, "_max")) end end for v in positivity_variables_cons