Skip to content

Commit

Permalink
Merge branch 'main' into sc/converters_coupling
Browse files Browse the repository at this point in the history
  • Loading branch information
SimonCan authored Oct 20, 2023
2 parents e730c9a + fc00ac1 commit 29f82f9
Show file tree
Hide file tree
Showing 9 changed files with 203 additions and 8 deletions.
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 @@ -235,7 +235,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
130 changes: 130 additions & 0 deletions src/callbacks_stage/subcell_bounds_check.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,130 @@
# 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)
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.
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
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,
solver.volume_integral, t,
iter + 1,
callback.output_directory,
save_errors)
end

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::BoundsCheckCallback, semi)
init_callback(callback, semi, semi.solver.volume_integral)
end

init_callback(callback::BoundsCheckCallback, semi, volume_integral::AbstractVolumeIntegral) = nothing

function init_callback(callback::BoundsCheckCallback, 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 v in limiter.positivity_variables_cons
print(f, ", " * string(variables[v]) * "_min")
end
end
println(f)
end

return nothing
end

function finalize_callback(callback::BoundsCheckCallback, semi)
finalize_callback(callback, semi, semi.solver.volume_integral)
end

finalize_callback(callback::BoundsCheckCallback, semi, volume_integral::AbstractVolumeIntegral) = nothing

function finalize_callback(callback::BoundsCheckCallback, 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(string(variables[v]) * ":\n- positivity: ",
idp_bounds_delta[Symbol(string(v), "_min")][2])
end
end
println(""^100 * "\n")

return nothing
end

include("subcell_bounds_check_2d.jl")
end # @muladd
49 changes: 49 additions & 0 deletions src/callbacks_stage/subcell_bounds_check_2d.jl
Original file line number Diff line number Diff line change
@@ -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)
(; positivity) = solver.volume_integral.limiter
(; variable_bounds) = limiter.cache.subcell_limiter_coefficients
(; idp_bounds_delta) = limiter.cache

if positivity
for v in limiter.positivity_variables_cons
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[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(string(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
end
end # @muladd
4 changes: 2 additions & 2 deletions 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 All @@ -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)

Expand Down
12 changes: 10 additions & 2 deletions src/solvers/dgsem_tree/subcell_limiters_2d.jl
Original file line number Diff line number Diff line change
Expand Up @@ -13,7 +13,15 @@ function create_cache(limiter::Type{SubcellLimiterIDP}, equations::AbstractEquat
nnodes(basis),
bound_keys)

return (; subcell_limiter_coefficients)
# 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)
end

return (; subcell_limiter_coefficients, idp_bounds_delta)
end

function (limiter::SubcellLimiterIDP)(u::AbstractArray{<:Any, 4}, semi, dg::DGSEM, t,
Expand Down Expand Up @@ -60,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]
Expand Down
7 changes: 6 additions & 1 deletion test/test_tree_2d_eulermulti.jl
Original file line number Diff line number Diff line change
Expand Up @@ -20,11 +20,16 @@ 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],
initial_refinement_level = 3,
tspan = (0.0, 0.001))
tspan = (0.0, 0.001),
output_directory="out")
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
Expand Down

0 comments on commit 29f82f9

Please sign in to comment.