Skip to content

Commit

Permalink
Add revised BoundsCheck for IDP Limiting
Browse files Browse the repository at this point in the history
  • Loading branch information
@bennibolm
bennibolm committed Oct 22, 2023
1 parent 3848564 commit 603d0f0
Show file tree
Hide file tree
Showing 3 changed files with 101 additions and 107 deletions.
67 changes: 43 additions & 24 deletions src/callbacks_stage/subcell_bounds_check.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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) &&
Expand All @@ -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

Expand All @@ -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
Expand All @@ -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)
Expand All @@ -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
Expand All @@ -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
Expand All @@ -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: ",
Expand All @@ -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)
Expand Down
137 changes: 56 additions & 81 deletions src/callbacks_stage/subcell_bounds_check_2d.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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)
Expand Down
4 changes: 2 additions & 2 deletions src/solvers/dgsem_tree/subcell_limiters_2d.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down

0 comments on commit 603d0f0

Please sign in to comment.