Merge branch 'main' into subcell-limiting-visualization

docs/src/callbacks.md

@@ -30,6 +30,11 @@ An example elixir using AMR can be found at [`examples/tree_2d_dgsem/elixir_adve
 The [`AnalysisCallback`](@ref) can be used to analyze the numerical solution, e.g. calculate
 errors or user-specified integrals, and print the results to the screen. The results can also be
 saved in a file. An example can be found at [`examples/tree_2d_dgsem/elixir_euler_vortex.jl`](https://github.com/trixi-framework/Trixi.jl/blob/main/examples/tree_2d_dgsem/elixir_euler_vortex.jl).
+Note that the errors (e.g. `L2 error` or `Linf error`) are computed with respect to the initial condition.
+The percentage of the simulation time refers to the ratio of the current time and the final time, i.e. it does
+not consider the maximal number of iterations. So the simulation could finish before 100% are reached.
+Note that, e.g., due to AMR or smaller time step sizes, the simulation can actually take longer than
+the percentage indicates.
 In [Performance metrics of the `AnalysisCallback`](@ref performance-metrics) you can find a detailed
 description of the different performance metrics the `AnalysisCallback` computes.
 

src/callbacks_step/alive.jl

@@ -86,9 +86,15 @@ function (alive_callback::AliveCallback)(integrator)
         println("─"^100)
         println()
     elseif mpi_isroot()
+        t = integrator.t
+        t_initial = first(integrator.sol.prob.tspan)
+        t_final = last(integrator.sol.prob.tspan)
+        sim_time_percentage = (t - t_initial) / (t_final - t_initial) * 100
         runtime_absolute = 1.0e-9 * (time_ns() - alive_callback.start_time)
-        @printf("#timesteps: %6d │ Δt: %.4e │ sim. time: %.4e │ run time: %.4e s\n",
-                integrator.stats.naccept, integrator.dt, integrator.t, runtime_absolute)
+        println(rpad(@sprintf("#timesteps: %6d │ Δt: %.4e │ sim. time: %.4e (%5.3f%%)",
+                              integrator.stats.naccept, integrator.dt, t,
+                              sim_time_percentage), 71) *
+                @sprintf("│ run time: %.4e s", runtime_absolute))
     end
 
     # avoid re-evaluating possible FSAL stages

src/callbacks_step/analysis.jl

@@ -232,6 +232,12 @@ function (analysis_callback::AnalysisCallback)(u_ode, du_ode, integrator, semi)
     @unpack dt, t = integrator
     iter = integrator.stats.naccept
 
+    # Compute the percentage of the simulation that is done
+    t = integrator.t
+    t_initial = first(integrator.sol.prob.tspan)
+    t_final = last(integrator.sol.prob.tspan)
+    sim_time_percentage = (t - t_initial) / (t_final - t_initial) * 100
+
     # Record performance measurements and compute performance index (PID)
     runtime_since_last_analysis = 1.0e-9 * (time_ns() -
                                    analysis_callback.start_time_last_analysis)
@@ -291,8 +297,8 @@ function (analysis_callback::AnalysisCallback)(u_ode, du_ode, integrator, semi)
                     "               " *
                     " └── GC time:    " *
                     @sprintf("%10.8e s (%5.3f%%)", gc_time_absolute, gc_time_percentage))
-        mpi_println(" sim. time:      " * @sprintf("%10.8e", t) *
-                    "               " *
+        mpi_println(rpad(" sim. time:      " *
+                         @sprintf("%10.8e (%5.3f%%)", t, sim_time_percentage), 46) *
                     " time/DOF/rhs!:  " * @sprintf("%10.8e s", runtime_relative))
         mpi_println("                 " * "              " *
                     "               " *

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,
@@ -1369,7 +1369,7 @@ nnodes(container::ContainerSubcellLimiterIDP2D) = size(container.alpha, 1)
 function Base.resize!(container::ContainerSubcellLimiterIDP2D, capacity)
     n_nodes = nnodes(container)
 
-    @unpack _alpha, _alpha1, _alpha2 = container
+    (; _alpha, _alpha1, _alpha2) = container
     resize!(_alpha, n_nodes * n_nodes * capacity)
     container.alpha = unsafe_wrap(Array, pointer(_alpha), (n_nodes, n_nodes, capacity))
     container.alpha .= convert(eltype(container.alpha), NaN)
@@ -1380,11 +1380,12 @@ function Base.resize!(container::ContainerSubcellLimiterIDP2D, capacity)
     container.alpha2 = unsafe_wrap(Array, pointer(_alpha2),
                                    (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))
+    (; _variable_bounds) = container
+    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/subcell_limiters.jl

@@ -52,9 +52,10 @@ function SubcellLimiterIDP(equations::AbstractEquations, basis;
                            positivity_variables_cons = [],
                            positivity_correction_factor = 0.1)
     positivity = (length(positivity_variables_cons) > 0)
-    number_bounds = length(positivity_variables_cons)
 
-    cache = create_cache(SubcellLimiterIDP, equations, basis, number_bounds)
+    bound_keys = Tuple(Symbol("$(i)_min") for i in positivity_variables_cons)
+
+    cache = create_cache(SubcellLimiterIDP, equations, basis, bound_keys)
 
     SubcellLimiterIDP{typeof(positivity_correction_factor), typeof(cache)}(positivity,
                                                                            positivity_variables_cons,

src/solvers/dgsem_tree/subcell_limiters_2d.jl

@@ -7,15 +7,13 @@
 
 # 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)
+                      basis::LobattoLegendreBasis, bound_keys)
     subcell_limiter_coefficients = Trixi.ContainerSubcellLimiterIDP2D{real(basis)
                                                                       }(0,
                                                                         nnodes(basis),
-                                                                        number_bounds)
+                                                                        bound_keys)
 
-    cache = (; subcell_limiter_coefficients)
-
-    return cache
+    return (; subcell_limiter_coefficients)
 end
 
 function (limiter::SubcellLimiterIDP)(u::AbstractArray{<:Any, 4}, semi, dg::DGSEM, t,
@@ -25,8 +23,7 @@ function (limiter::SubcellLimiterIDP)(u::AbstractArray{<:Any, 4}, semi, dg::DGSE
     alpha .= zero(eltype(alpha))
 
     if limiter.positivity
-        @trixi_timeit timer() "positivity" idp_positivity!(alpha, limiter, u, dt,
-                                                           semi)
+        @trixi_timeit timer() "positivity" idp_positivity!(alpha, limiter, u, dt, semi)
     end
 
     # Calculate alpha1 and alpha2
@@ -49,22 +46,21 @@ end
 
 @inline function idp_positivity!(alpha, limiter, u, dt, semi)
     # Conservative variables
-    for (index, variable) in enumerate(limiter.positivity_variables_cons)
-        idp_positivity!(alpha, limiter, u, dt, semi, variable, index)
+    for variable in limiter.positivity_variables_cons
+        idp_positivity!(alpha, limiter, u, dt, semi, variable)
     end
 
     return nothing
 end
 
-@inline function idp_positivity!(alpha, limiter, u, dt, semi, variable, index)
+@inline function idp_positivity!(alpha, limiter, u, dt, semi, variable)
     mesh, equations, dg, cache = mesh_equations_solver_cache(semi)
-    @unpack antidiffusive_flux1, antidiffusive_flux2 = cache.antidiffusive_fluxes
-    @unpack inverse_weights = dg.basis
-    @unpack positivity_correction_factor = limiter
-
-    @unpack variable_bounds = limiter.cache.subcell_limiter_coefficients
+    (; antidiffusive_flux1, antidiffusive_flux2) = cache.antidiffusive_fluxes
+    (; inverse_weights) = dg.basis
+    (; positivity_correction_factor) = limiter
 
-    var_min = variable_bounds[index]
+    (; variable_bounds) = limiter.cache.subcell_limiter_coefficients
+    var_min = variable_bounds[Symbol("$(variable)_min")]
 
     @threaded for element in eachelement(dg, semi.cache)
         inverse_jacobian = cache.elements.inverse_jacobian[element]