Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add implementation of subcell bounds check #1672

Merged
merged 20 commits into from
Oct 20, 2023
Merged
Show file tree
Hide file tree
Changes from 6 commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -132,7 +132,9 @@ callbacks = CallbackSet(summary_callback,
###############################################################################
# run the simulation

stage_callbacks = (SubcellLimiterIDPCorrection(),)
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
Expand Down
2 changes: 1 addition & 1 deletion src/Trixi.jl
Original file line number Diff line number Diff line change
Expand Up @@ -231,7 +231,7 @@ export DG,
SurfaceIntegralUpwind,
MortarL2

export VolumeIntegralSubcellLimiting,
export VolumeIntegralSubcellLimiting, BoundsCheckCallback,
SubcellLimiterIDP, SubcellLimiterIDPCorrection

export nelements, nnodes, nvariables,
Expand Down
1 change: 1 addition & 0 deletions src/callbacks_stage/callbacks_stage.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
129 changes: 129 additions & 0 deletions src/callbacks_stage/subcell_bounds_check.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,129 @@
# 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.
bennibolm marked this conversation as resolved.
Show resolved Hide resolved

!!! 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
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)) &&
(iter % callback.interval == 0 || integrator.finalstep)
@trixi_timeit timer() "check_bounds" check_bounds(u, mesh, equations, solver, cache,
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
bennibolm marked this conversation as resolved.
Show resolved Hide resolved

function check_bounds(u, mesh, equations, solver, cache,
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)
check_bounds(u, mesh, equations, solver, cache, volume_integral.limiter, t, iter,
output_directory, save_errors)
end

function init_callback(callback, semi)
init_callback(callback, semi, semi.solver.volume_integral)
end

init_callback(callback, semi, volume_integral::AbstractVolumeIntegral) = nothing
bennibolm marked this conversation as resolved.
Show resolved Hide resolved

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
48 changes: 48 additions & 0 deletions src/callbacks_stage/subcell_bounds_check_2d.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,48 @@
# 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)
(; positivity) = 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 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
bennibolm marked this conversation as resolved.
Show resolved Hide resolved

return nothing
end
end # @muladd
2 changes: 1 addition & 1 deletion src/solvers/dgsem_tree/subcell_limiters.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down
7 changes: 6 additions & 1 deletion src/solvers/dgsem_tree/subcell_limiters_2d.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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)}()
bennibolm marked this conversation as resolved.
Show resolved Hide resolved
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,
Expand Down
6 changes: 5 additions & 1 deletion test/test_tree_2d_eulermulti.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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] == "15, 0.0009595796113231045, 0.0, 0.0"
bennibolm marked this conversation as resolved.
Show resolved Hide resolved
end

@trixi_testset "elixir_eulermulti_ec.jl" begin
Expand Down
Loading