From 2e251852d79c98e2328d0b97e79ad6cf9052f4a5 Mon Sep 17 00:00:00 2001 From: bennibolm Date: Thu, 28 Sep 2023 17:53:38 +0200 Subject: [PATCH] Implement variable_bounds as Dict with Symbols --- src/callbacks_stage/bounds_check.jl | 28 +++--- src/callbacks_stage/bounds_check_2d.jl | 51 +++++----- src/solvers/dgsem_tree/containers_2d.jl | 27 +++--- .../dgsem_tree/dg_2d_subcell_limiters.jl | 11 +-- src/solvers/dgsem_tree/subcell_limiters.jl | 27 ++++-- src/solvers/dgsem_tree/subcell_limiters_2d.jl | 93 +++++++------------ 6 files changed, 106 insertions(+), 131 deletions(-) diff --git a/src/callbacks_stage/bounds_check.jl b/src/callbacks_stage/bounds_check.jl index 1d56d725c8d..600c11d9252 100644 --- a/src/callbacks_stage/bounds_check.jl +++ b/src/callbacks_stage/bounds_check.jl @@ -145,34 +145,30 @@ end println("─"^100) println("Maximum deviation from bounds:") println("─"^100) - counter = 1 if local_minmax - for index in limiter.local_minmax_variables_cons - println("$(variables[index]):") - println("-lower bound: ", idp_bounds_delta[counter]) - println("-upper bound: ", idp_bounds_delta[counter + 1]) - counter += 2 + for v in limiter.local_minmax_variables_cons + println("$(variables[v]):") + println("-lower bound: ", idp_bounds_delta[Symbol("$(v)_min")]) + println("-upper bound: ", idp_bounds_delta[Symbol("$(v)_max")]) end end if spec_entropy - println("spec. entropy:\n- lower bound: ", idp_bounds_delta[counter]) - counter += 1 + println("spec. entropy:\n- lower bound: ", idp_bounds_delta[:spec_entropy_min]) end if math_entropy - println("math. entropy:\n- upper bound: ", idp_bounds_delta[counter]) - counter += 1 + println("math. entropy:\n- upper bound: ", idp_bounds_delta[:math_entropy_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 - println("$(variables[index]):\n- positivity: ", idp_bounds_delta[counter]) - counter += 1 + println("$(variables[v]):\n- positivity: ", + idp_bounds_delta[Symbol("$(v)_min")]) end for variable in limiter.positivity_variables_nonlinear - println("$(variable):\n- positivity: ", idp_bounds_delta[counter]) - counter += 1 + println("$(variable):\n- positivity: ", + idp_bounds_delta[Symbol("$(variable)_min")]) end end println("─"^100 * "\n") diff --git a/src/callbacks_stage/bounds_check_2d.jl b/src/callbacks_stage/bounds_check_2d.jl index 7b2740ff180..f66f8b8b3aa 100644 --- a/src/callbacks_stage/bounds_check_2d.jl +++ b/src/callbacks_stage/bounds_check_2d.jl @@ -13,29 +13,31 @@ @unpack idp_bounds_delta = limiter.cache save_errors_ = save_errors && (iter % interval == 0) - counter = 1 if save_errors_ open("$output_directory/deviations.txt", "a") do f print(f, iter, ", ", time) end end if local_minmax - for index in limiter.local_minmax_variables_cons + 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)) for element in eachelement(solver, cache), j in eachnode(solver), i in eachnode(solver) deviation_min = max(deviation_min, - variable_bounds[counter][i, j, element] - - u[index, i, j, element]) + variable_bounds[key_min][i, j, element] - + u[v, i, j, element]) deviation_max = max(deviation_max, - u[index, i, j, element] - - variable_bounds[counter + 1][i, j, element]) + u[v, i, j, element] - + variable_bounds[key_max][i, j, element]) end - idp_bounds_delta[counter] = max(idp_bounds_delta[counter], deviation_min) - idp_bounds_delta[counter + 1] = max(idp_bounds_delta[counter + 1], - deviation_max) + 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 @@ -43,10 +45,10 @@ print(f, ", ", deviation_min_, ", ", deviation_max_) end end - counter += 2 end end if spec_entropy + key = :spec_entropy_min deviation_min = zero(eltype(u)) for element in eachelement(solver, cache), j in eachnode(solver), i in eachnode(solver) @@ -54,18 +56,18 @@ s = entropy_spec(get_node_vars(u, equations, solver, i, j, element), equations) deviation_min = max(deviation_min, - variable_bounds[counter][i, j, element] - s) + variable_bounds[key][i, j, element] - s) end - idp_bounds_delta[counter] = max(idp_bounds_delta[counter], deviation_min) + 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 - counter += 1 end if math_entropy + key = :math_entropy_max deviation_max = zero(eltype(u)) for element in eachelement(solver, cache), j in eachnode(solver), i in eachnode(solver) @@ -73,40 +75,40 @@ s = entropy_math(get_node_vars(u, equations, solver, i, j, element), equations) deviation_max = max(deviation_max, - s - variable_bounds[counter][i, j, element]) + s - variable_bounds[key][i, j, element]) end - idp_bounds_delta[counter] = max(idp_bounds_delta[counter], deviation_max) + 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 end - counter += 1 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 + 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[index, i, j, element] + var = u[v, i, j, element] deviation_min = max(deviation_min, - variable_bounds[counter][i, j, element] - var) + variable_bounds[key][i, j, element] - var) end - idp_bounds_delta[counter] = max(idp_bounds_delta[counter], deviation_min) + 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 - counter += 1 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) @@ -114,16 +116,15 @@ var = variable(get_node_vars(u, equations, solver, i, j, element), equations) deviation_min = max(deviation_min, - variable_bounds[counter][i, j, element] - var) + variable_bounds[key][i, j, element] - var) end - idp_bounds_delta[counter] = max(idp_bounds_delta[counter], deviation_min) + 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 - counter += 1 end end if save_errors_ diff --git a/src/solvers/dgsem_tree/containers_2d.jl b/src/solvers/dgsem_tree/containers_2d.jl index b1bd509fc7d..a6d9f57758b 100644 --- a/src/solvers/dgsem_tree/containers_2d.jl +++ b/src/solvers/dgsem_tree/containers_2d.jl @@ -1325,16 +1325,16 @@ mutable struct ContainerSubcellLimiterIDP2D{uEltype <: Real} alpha::Array{uEltype, 3} # [i, j, element] alpha1::Array{uEltype, 3} alpha2::Array{uEltype, 3} - variable_bounds::Vector{Array{uEltype, 3}} + variable_bounds::Dict{Symbol, Array{uEltype, 3}} # internal `resize!`able storage _alpha::Vector{uEltype} _alpha1::Vector{uEltype} _alpha2::Vector{uEltype} - _variable_bounds::Vector{Vector{uEltype}} + _variable_bounds::Dict{Symbol, Vector{uEltype}} end function ContainerSubcellLimiterIDP2D{uEltype}(capacity::Integer, n_nodes, - length) where {uEltype <: Real} + bound_keys) where {uEltype <: Real} nan_uEltype = convert(uEltype, NaN) # Initialize fields with defaults @@ -1345,12 +1345,12 @@ function ContainerSubcellLimiterIDP2D{uEltype}(capacity::Integer, n_nodes, _alpha2 = fill(nan_uEltype, n_nodes * (n_nodes + 1) * capacity) alpha2 = unsafe_wrap(Array, pointer(_alpha2), (n_nodes, n_nodes + 1, capacity)) - _variable_bounds = Vector{Vector{uEltype}}(undef, length) - variable_bounds = Vector{Array{uEltype, 3}}(undef, length) - for i in 1:length - _variable_bounds[i] = fill(nan_uEltype, n_nodes * n_nodes * capacity) - variable_bounds[i] = unsafe_wrap(Array, pointer(_variable_bounds[i]), - (n_nodes, n_nodes, capacity)) + _variable_bounds = Dict{Symbol, Vector{uEltype}}() + variable_bounds = Dict{Symbol, Array{uEltype, 3}}() + for key in bound_keys + _variable_bounds[key] = fill(nan_uEltype, n_nodes * n_nodes * capacity) + variable_bounds[key] = unsafe_wrap(Array, pointer(_variable_bounds[key]), + (n_nodes, n_nodes, capacity)) end return ContainerSubcellLimiterIDP2D{uEltype}(alpha, alpha1, alpha2, @@ -1380,10 +1380,11 @@ function Base.resize!(container::ContainerSubcellLimiterIDP2D, capacity) (n_nodes, n_nodes + 1, capacity)) @unpack _variable_bounds = container - for i in 1:length(_variable_bounds) - resize!(_variable_bounds[i], n_nodes * n_nodes * capacity) - container.variable_bounds[i] = unsafe_wrap(Array, pointer(_variable_bounds[i]), - (n_nodes, n_nodes, capacity)) + for (key, _) in _variable_bounds + resize!(_variable_bounds[key], n_nodes * n_nodes * capacity) + container.variable_bounds[key] = unsafe_wrap(Array, + pointer(_variable_bounds[key]), + (n_nodes, n_nodes, capacity)) end return nothing diff --git a/src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl b/src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl index 0f0fb1a4b71..e165308967f 100644 --- a/src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl +++ b/src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl @@ -534,12 +534,11 @@ end @unpack variable_bounds = limiter.cache.subcell_limiter_coefficients @unpack bar_states1, bar_states2 = limiter.cache.container_bar_states - counter = 1 # state variables if limiter.local_minmax for index in limiter.local_minmax_variables_cons - var_min = variable_bounds[counter] - var_max = variable_bounds[counter + 1] + var_min = variable_bounds[Symbol("$(index)_min")] + var_max = variable_bounds[Symbol("$(index)_max")] @threaded for element in eachelement(dg, cache) var_min[:, :, element] .= typemax(eltype(var_min)) var_max[:, :, element] .= typemin(eltype(var_max)) @@ -571,12 +570,11 @@ end bar_states2[index, i, j + 1, element]) end end - counter += 2 end end # Specific Entropy if limiter.spec_entropy - s_min = variable_bounds[counter] + s_min = variable_bounds[:spec_entropy_min] @threaded for element in eachelement(dg, cache) s_min[:, :, element] .= typemax(eltype(s_min)) for j in eachnode(dg), i in eachnode(dg) @@ -602,11 +600,10 @@ end s_min[i, j, element] = min(s_min[i, j, element], s) end end - counter += 1 end # Mathematical entropy if limiter.math_entropy - s_max = variable_bounds[counter] + s_max = variable_bounds[:math_entropy_max] @threaded for element in eachelement(dg, cache) s_max[:, :, element] .= typemin(eltype(s_max)) for j in eachnode(dg), i in eachnode(dg) diff --git a/src/solvers/dgsem_tree/subcell_limiters.jl b/src/solvers/dgsem_tree/subcell_limiters.jl index 8dd17c2d6a9..86d88c827f8 100644 --- a/src/solvers/dgsem_tree/subcell_limiters.jl +++ b/src/solvers/dgsem_tree/subcell_limiters.jl @@ -108,17 +108,28 @@ function SubcellLimiterIDP(equations::AbstractEquations, basis; error("Only one of the two can be selected: math_entropy/spec_entropy") end - number_bounds = 2 * length(local_minmax_variables_cons) + - length(positivity_variables_nonlinear) + - spec_entropy + math_entropy - - for index in positivity_variables_cons - if !(index in local_minmax_variables_cons) - number_bounds += 1 + bound_keys = () + if local_minmax + for i in local_minmax_variables_cons + bound_keys = (bound_keys..., Symbol("$(i)_min"), Symbol("$(i)_max")) + end + end + if spec_entropy + bound_keys = (bound_keys..., :spec_entropy_min) + end + if math_entropy + bound_keys = (bound_keys..., :math_entropy_max) + end + for i in positivity_variables_cons + if !(i in local_minmax_variables_cons) + bound_keys = (bound_keys..., Symbol("$(i)_min")) end end + for variable in positivity_variables_nonlinear + bound_keys = (bound_keys..., Symbol("$(variable)_min")) + end - cache = create_cache(SubcellLimiterIDP, equations, basis, number_bounds, bar_states) + cache = create_cache(SubcellLimiterIDP, equations, basis, bound_keys, bar_states) if smoothness_indicator IndicatorHG = IndicatorHennemannGassner(equations, basis, alpha_max = 1.0, diff --git a/src/solvers/dgsem_tree/subcell_limiters_2d.jl b/src/solvers/dgsem_tree/subcell_limiters_2d.jl index 7db04ca83ff..a8fcc42c061 100644 --- a/src/solvers/dgsem_tree/subcell_limiters_2d.jl +++ b/src/solvers/dgsem_tree/subcell_limiters_2d.jl @@ -7,11 +7,11 @@ # this method is used when the limiter is constructed as for shock-capturing volume integrals function create_cache(limiter::Type{SubcellLimiterIDP}, equations::AbstractEquations{2}, - basis::LobattoLegendreBasis, number_bounds, bar_states) + basis::LobattoLegendreBasis, bound_keys, bar_states) subcell_limiter_coefficients = Trixi.ContainerSubcellLimiterIDP2D{real(basis) }(0, nnodes(basis), - number_bounds) + bound_keys) cache = (;) if bar_states @@ -21,7 +21,10 @@ function create_cache(limiter::Type{SubcellLimiterIDP}, equations::AbstractEquat cache = (; cache..., container_bar_states) end - idp_bounds_delta = zeros(real(basis), number_bounds) + idp_bounds_delta = Dict{Symbol, real(basis)}() + for key in bound_keys + idp_bounds_delta[key] = zero(real(basis)) + end return (; cache..., subcell_limiter_coefficients, idp_bounds_delta) end @@ -318,20 +321,18 @@ end end @inline function idp_local_minmax!(alpha, limiter, u, t, dt, semi, elements) - for (index, variable) in enumerate(limiter.local_minmax_variables_cons) - idp_local_minmax!(alpha, limiter, u, t, dt, semi, elements, variable, index) + for variable in limiter.local_minmax_variables_cons + idp_local_minmax!(alpha, limiter, u, t, dt, semi, elements, variable) end return nothing end -@inline function idp_local_minmax!(alpha, limiter, u, t, dt, semi, elements, variable, - index) +@inline function idp_local_minmax!(alpha, limiter, u, t, dt, semi, elements, variable) mesh, _, dg, cache = mesh_equations_solver_cache(semi) - @unpack variable_bounds = limiter.cache.subcell_limiter_coefficients - - var_min = variable_bounds[2 * (index - 1) + 1] - var_max = variable_bounds[2 * (index - 1) + 2] + (; variable_bounds) = limiter.cache.subcell_limiter_coefficients + var_min = variable_bounds[Symbol("$(variable)_min")] + var_max = variable_bounds[Symbol("$(variable)_max")] if !limiter.bar_states calc_bounds_2sided!(var_min, var_max, variable, u, t, semi) end @@ -396,9 +397,8 @@ end @inline function idp_spec_entropy!(alpha, limiter, u, t, dt, semi, elements) mesh, equations, dg, cache = mesh_equations_solver_cache(semi) - @unpack variable_bounds = limiter.cache.subcell_limiter_coefficients - - s_min = variable_bounds[2 * length(limiter.local_minmax_variables_cons) + 1] + (; variable_bounds) = limiter.cache.subcell_limiter_coefficients + s_min = variable_bounds[:spec_entropy_min] if !limiter.bar_states calc_bounds_1sided!(s_min, min, typemax, entropy_spec, u, t, semi) end @@ -427,10 +427,8 @@ end @inline function idp_math_entropy!(alpha, limiter, u, t, dt, semi, elements) mesh, equations, dg, cache = mesh_equations_solver_cache(semi) - @unpack spec_entropy = limiter - @unpack variable_bounds = limiter.cache.subcell_limiter_coefficients - - s_max = variable_bounds[2 * length(limiter.local_minmax_variables_cons) + spec_entropy + 1] + (; variable_bounds) = limiter.cache.subcell_limiter_coefficients + s_max = variable_bounds[:math_entropy_max] if !limiter.bar_states calc_bounds_1sided!(s_max, max, typemin, entropy_math, u, t, semi) end @@ -459,48 +457,26 @@ end @inline function idp_positivity!(alpha, limiter, u, dt, semi, elements) # Conservative variables - for (index, variable) in enumerate(limiter.positivity_variables_cons) - idp_positivity!(alpha, limiter, u, dt, semi, elements, variable, index) + for variable in limiter.positivity_variables_cons + idp_positivity!(alpha, limiter, u, dt, semi, elements, variable) end # Nonlinear variables - for (index, variable) in enumerate(limiter.positivity_variables_nonlinear) - idp_positivity_newton!(alpha, limiter, u, dt, semi, elements, variable, index) + for variable in limiter.positivity_variables_nonlinear + idp_positivity_nonlinear!(alpha, limiter, u, dt, semi, elements, variable) end return nothing end -@inline function idp_positivity!(alpha, limiter, u, dt, semi, elements, variable, - index) +@inline function idp_positivity!(alpha, limiter, u, dt, semi, elements, variable) mesh, equations, dg, cache = mesh_equations_solver_cache(semi) @unpack antidiffusive_flux1, antidiffusive_flux2 = cache.antidiffusive_fluxes @unpack inverse_weights = dg.basis - @unpack local_minmax, spec_entropy, math_entropy, positivity_correction_factor = limiter - - @unpack variable_bounds = limiter.cache.subcell_limiter_coefficients - - counter = 2 * length(limiter.local_minmax_variables_cons) + spec_entropy + - math_entropy - if local_minmax - if variable in limiter.local_minmax_variables_cons - for (index_, variable_) in enumerate(limiter.local_minmax_variables_cons) - if variable == variable_ - var_min = variable_bounds[2 * (index_ - 1) + 1] - break - end - end - else - for variable_ in limiter.positivity_variables_cons[1:index] - if !(variable_ in limiter.local_minmax_variables_cons) - counter += 1 - end - end - var_min = variable_bounds[counter] - end - else - var_min = variable_bounds[counter + index] - end + @unpack local_minmax, positivity_correction_factor = limiter + + (; variable_bounds) = limiter.cache.subcell_limiter_coefficients + var_min = variable_bounds[Symbol("$(variable)_min")] @threaded for element in elements if mesh isa TreeMesh @@ -558,20 +534,13 @@ end return nothing end -@inline function idp_positivity_newton!(alpha, limiter, u, dt, semi, elements, - variable, index) +@inline function idp_positivity_nonlinear!(alpha, limiter, u, dt, semi, elements, + variable) mesh, equations, dg, cache = mesh_equations_solver_cache(semi) - @unpack spec_entropy, math_entropy, positivity_correction_factor, positivity_variables_cons = limiter - @unpack variable_bounds = limiter.cache.subcell_limiter_coefficients - - index_ = 2 * length(limiter.local_minmax_variables_cons) + spec_entropy + - math_entropy + index - for variable_ in limiter.positivity_variables_cons - if !(variable_ in limiter.local_minmax_variables_cons) - index_ += 1 - end - end - var_min = variable_bounds[index_] + (; positivity_correction_factor) = limiter + + (; variable_bounds) = limiter.cache.subcell_limiter_coefficients + var_min = variable_bounds[Symbol("$(variable)_min")] @threaded for element in elements for j in eachnode(dg), i in eachnode(dg)