Merge pull request #114 from bennibolm/subcell-limiting-dictionary

Implement variable_bounds as Dict

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")

src/callbacks_stage/bounds_check_2d.jl

@@ -13,117 +13,118 @@
     @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
                 open("$output_directory/deviations.txt", "a") do f
                     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)
 
             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)
 
             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)
 
                 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_

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

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)

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,