diff --git a/src/callbacks_stage/subcell_bounds_check.jl b/src/callbacks_stage/subcell_bounds_check.jl index 34507dced50..9a033e72367 100644 --- a/src/callbacks_stage/subcell_bounds_check.jl +++ b/src/callbacks_stage/subcell_bounds_check.jl @@ -8,9 +8,19 @@ """ BoundsCheckCallback(; output_directory="out", save_errors=false, interval=1) -Bounds checking routine for [`SubcellLimiterIDP`](@ref) and [`SubcellLimiterMCL`](@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) and [`SubcellLimiterMCL`](@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. +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". +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 + [`SubcellLimiterIDPCorrection`](@ref). So, to check the final solution, this bounds check + callback must be called after the correction stage. """ struct BoundsCheckCallback output_directory::String @@ -25,7 +35,7 @@ end function (callback::BoundsCheckCallback)(u_ode, integrator, stage) mesh, equations, solver, cache = mesh_equations_solver_cache(integrator.p) - @unpack t, iter, alg = integrator + (; t, iter, alg) = integrator u = wrap_array(u_ode, mesh, equations, solver, cache) save_errors = callback.save_errors && (callback.interval > 0) && @@ -51,13 +61,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 @@ -66,16 +77,17 @@ function init_callback(callback::BoundsCheckCallback, semi, limiter::SubcellLimi return nothing end - @unpack local_minmax, positivity, spec_entropy, math_entropy = limiter - @unpack output_directory = callback + (; local_minmax, positivity, spec_entropy, math_entropy) = 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 index in limiter.local_minmax_variables_cons - print(f, ", $(variables[index])_min, $(variables[index])_max") + for v in limiter.local_minmax_variables_cons + variable_string = string(variables[v]) + print(f, ", " * variable_string * "_min, " * variable_string * "_max") end end if spec_entropy @@ -85,14 +97,14 @@ function init_callback(callback::BoundsCheckCallback, semi, limiter::SubcellLimi print(f, ", mathEntr_max") end if positivity - for index in limiter.positivity_variables_cons - if index in limiter.local_minmax_variables_cons + for v in limiter.positivity_variables_cons + if v in limiter.local_minmax_variables_cons continue end - print(f, ", $(variables[index])_min") + print(f, ", " * string(variables[v]) * "_min") end for variable in limiter.positivity_variables_nonlinear - print(f, ", $(variable)_min") + print(f, ", " * string(variable) * "_min") end end println(f) @@ -106,6 +118,8 @@ function init_callback(callback::BoundsCheckCallback, semi, limiter::SubcellLimi return nothing end + # TODO: Revise Bounds Check for MCL + @unpack output_directory = callback mkpath(output_directory) open("$output_directory/deviations.txt", "a") do f @@ -121,13 +135,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 @@ -143,24 +157,27 @@ end 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)_min")]) - println("-upper bound: ", idp_bounds_delta[Symbol("$(v)_max")]) + 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 spec_entropy - println("spec. entropy:\n- lower bound: ", idp_bounds_delta[:spec_entropy_min]) + println("spec. entropy:\n- lower bound: ", + idp_bounds_delta[:spec_entropy_min][2]) end if math_entropy - println("math. entropy:\n- upper bound: ", idp_bounds_delta[:math_entropy_max]) + println("math. entropy:\n- upper bound: ", + idp_bounds_delta[:math_entropy_max][2]) end if positivity for v in limiter.positivity_variables_cons if v in limiter.local_minmax_variables_cons continue end - println("$(variables[v]):\n- positivity: ", - idp_bounds_delta[Symbol("$(v)_min")]) + println(string(variables[v]) * ":\n- positivity: ", + idp_bounds_delta[Symbol(string(v), "_min")][2]) end for variable in limiter.positivity_variables_nonlinear println("$(variable):\n- positivity: ", @@ -176,6 +193,8 @@ end limiter::SubcellLimiterMCL) @unpack idp_bounds_delta = limiter.cache + # TODO: Revise bounds check for MCL + println("─"^100) println("Maximum deviation from bounds:") println("─"^100) diff --git a/src/callbacks_stage/subcell_bounds_check_2d.jl b/src/callbacks_stage/subcell_bounds_check_2d.jl index 8334c108cf9..fe5be3d02b2 100644 --- a/src/callbacks_stage/subcell_bounds_check_2d.jl +++ b/src/callbacks_stage/subcell_bounds_check_2d.jl @@ -8,128 +8,101 @@ @inline function check_bounds(u, mesh::AbstractMesh{2}, equations, solver, cache, limiter::SubcellLimiterIDP, time, iter, output_directory, save_errors) - @unpack local_minmax, positivity, spec_entropy, math_entropy = solver.volume_integral.limiter - @unpack variable_bounds = limiter.cache.subcell_limiter_coefficients - @unpack idp_bounds_delta = limiter.cache + (; local_minmax, positivity, spec_entropy, math_entropy) = 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 local_minmax for v in limiter.local_minmax_variables_cons - key_min = Symbol("$(v)_min") - key_max = Symbol("$(v)_max") - deviation_min = zero(eltype(u)) - deviation_max = zero(eltype(u)) + 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) - deviation_min = max(deviation_min, - variable_bounds[key_min][i, j, element] - - u[v, i, j, element]) - deviation_max = max(deviation_max, - u[v, i, j, element] - - variable_bounds[key_max][i, j, element]) - end - idp_bounds_delta[key_min] = max(idp_bounds_delta[key_min], - deviation_min) - idp_bounds_delta[key_max] = max(idp_bounds_delta[key_max], - deviation_max) - if save_errors - deviation_min_ = deviation_min - deviation_max_ = deviation_max - open("$output_directory/deviations.txt", "a") do f - print(f, ", ", deviation_min_, ", ", deviation_max_) - end + 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 spec_entropy key = :spec_entropy_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) s = entropy_spec(get_node_vars(u, equations, solver, i, j, element), equations) - deviation_min = max(deviation_min, - variable_bounds[key][i, j, element] - s) - 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] - s) end + deviation[2] = max(deviation[2], deviation[1]) end if math_entropy key = :math_entropy_max - deviation_max = zero(eltype(u)) + deviation = idp_bounds_delta[key] for element in eachelement(solver, cache), j in eachnode(solver), i in eachnode(solver) s = entropy_math(get_node_vars(u, equations, solver, i, j, element), equations) - deviation_max = max(deviation_max, - s - variable_bounds[key][i, j, element]) - end - idp_bounds_delta[key] = max(idp_bounds_delta[key], deviation_max) - if save_errors - deviation_max_ = deviation_max - open("$output_directory/deviations.txt", "a") do f - print(f, ", ", deviation_max_) - end + deviation[1] = max(deviation[1], s - variable_bounds[key][i, j, element]) end + deviation[2] = max(deviation[2], deviation[1]) end if positivity for v in limiter.positivity_variables_cons - if v in limiter.local_minmax_variables_cons - continue - end - key = Symbol("$(v)_min") - deviation_min = zero(eltype(u)) + key = Symbol(string(v), "_min") + 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) + deviation[1] = max(deviation[1], + 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_) + 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 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") + print(f, ", ", idp_bounds_delta[key_min][1], + idp_bounds_delta[key_max][1]) end end - end - for variable in limiter.positivity_variables_nonlinear - key = Symbol("$(variable)_min") - deviation_min = zero(eltype(u)) - for element in eachelement(solver, cache), j in eachnode(solver), - i in eachnode(solver) - - var = variable(get_node_vars(u, equations, solver, i, j, element), - equations) - deviation_min = max(deviation_min, - variable_bounds[key][i, j, element] - var) + if math_entropy + key = :math_entropy_max + print(f, ", ", idp_bounds_delta[key][1]) + end + if math_entropy + key = :math_entropy_max + print(f, ", ", idp_bounds_delta[key][1]) 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_) + if positivity + for v in limiter.positivity_variables_cons + key = Symbol(string(v), "_min") + print(f, ", ", idp_bounds_delta[key][1]) end end - end - end - if save_errors - open("$output_directory/deviations.txt", "a") do f 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 @@ -143,6 +116,8 @@ end @unpack idp_bounds_delta = limiter.cache @unpack antidiffusive_flux1, antidiffusive_flux2 = cache.antidiffusive_fluxes + # TODO: Revise Bounds Check for MCL + n_vars = nvariables(equations) deviation_min = zeros(eltype(u), n_vars + limiter.PressurePositivityLimiterKuzmin) diff --git a/src/solvers/dgsem_tree/subcell_limiters_2d.jl b/src/solvers/dgsem_tree/subcell_limiters_2d.jl index 83459c108b4..19a744f38d5 100644 --- a/src/solvers/dgsem_tree/subcell_limiters_2d.jl +++ b/src/solvers/dgsem_tree/subcell_limiters_2d.jl @@ -24,9 +24,9 @@ function create_cache(limiter::Type{SubcellLimiterIDP}, equations::AbstractEquat # 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, 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 (; cache..., subcell_limiter_coefficients, idp_bounds_delta)