Node-level visualization support for coefficients of Subcell limiting (…

…#1611)

* Add node-level visualization for IDP limiting

* Fix renaming error

* Init alphas with 0 after resizing

* Remove note

* Relocate overwriting of alphas

* Implement suggestions

* Add comment about split into nominator/denominator

* Set uninitialized alphas to NaN

src/callbacks_step/save_solution.jl

@@ -222,21 +222,28 @@ end
         end
     end
 
+    node_variables = Dict{Symbol, Any}()
+    @trixi_timeit timer() "get node variables" get_node_variables!(node_variables,
+                                                                   semi)
+
     @trixi_timeit timer() "save solution" save_solution_file(u_ode, t, dt, iter, semi,
                                                              solution_callback,
                                                              element_variables,
+                                                             node_variables,
                                                              system = system)
 end
 
 @inline function save_solution_file(u_ode, t, dt, iter,
                                     semi::AbstractSemidiscretization, solution_callback,
-                                    element_variables = Dict{Symbol, Any}();
+                                    element_variables = Dict{Symbol, Any}(),
+                                    node_variables = Dict{Symbol, Any}();
                                     system = "")
     mesh, equations, solver, cache = mesh_equations_solver_cache(semi)
     u = wrap_array_native(u_ode, mesh, equations, solver, cache)
     save_solution_file(u, t, dt, iter, mesh, equations, solver, cache,
                        solution_callback,
-                       element_variables; system = system)
+                       element_variables,
+                       node_variables; system = system)
 end
 
 # TODO: Taal refactor, move save_mesh_file?

src/callbacks_step/save_solution_dg.jl

@@ -10,7 +10,9 @@ function save_solution_file(u, time, dt, timestep,
                                         UnstructuredMesh2D, SerialP4estMesh,
                                         SerialT8codeMesh},
                             equations, dg::DG, cache,
-                            solution_callback, element_variables = Dict{Symbol, Any}();
+                            solution_callback,
+                            element_variables = Dict{Symbol, Any}(),
+                            node_variables = Dict{Symbol, Any}();
                             system = "")
     @unpack output_directory, solution_variables = solution_callback
 
@@ -73,6 +75,16 @@ function save_solution_file(u, time, dt, timestep,
             var = file["element_variables_$v"]
             attributes(var)["name"] = string(key)
         end
+
+        # Store node variables
+        for (v, (key, node_variable)) in enumerate(node_variables)
+            # Add to file
+            file["node_variables_$v"] = node_variable
+
+            # Add variable name as attribute
+            var = file["node_variables_$v"]
+            attributes(var)["name"] = string(key)
+        end
     end
 
     return filename
@@ -81,7 +93,9 @@ end
 function save_solution_file(u, time, dt, timestep,
                             mesh::Union{ParallelTreeMesh, ParallelP4estMesh}, equations,
                             dg::DG, cache,
-                            solution_callback, element_variables = Dict{Symbol, Any}();
+                            solution_callback,
+                            element_variables = Dict{Symbol, Any}(),
+                            node_variables = Dict{Symbol, Any}();
                             system = "")
     @unpack output_directory, solution_variables = solution_callback
 

src/semidiscretization/semidiscretization.jl

@@ -335,6 +335,10 @@ function get_element_variables!(element_variables, u_ode,
     get_element_variables!(element_variables, u, mesh_equations_solver_cache(semi)...)
 end
 
+function get_node_variables!(node_variables, semi::AbstractSemidiscretization)
+    get_node_variables!(node_variables, mesh_equations_solver_cache(semi)...)
+end
+
 # To implement AMR and use OrdinaryDiffEq.jl etc., we have to be a bit creative.
 # Since the caches of the SciML ecosystem are immutable structs, we cannot simply
 # change the underlying arrays therein. Hence, to support changing the number of

src/solvers/dg.jl

@@ -12,6 +12,11 @@ function get_element_variables!(element_variables, u, mesh, equations,
     nothing
 end
 
+function get_node_variables!(node_variables, mesh, equations,
+                             volume_integral::AbstractVolumeIntegral, dg, cache)
+    nothing
+end
+
 """
     VolumeIntegralStrongForm()
 
@@ -214,6 +219,18 @@ function Base.show(io::IO, mime::MIME"text/plain",
     end
 end
 
+function get_node_variables!(node_variables, mesh, equations,
+                             volume_integral::VolumeIntegralSubcellLimiting, dg, cache)
+    # While for the element-wise limiting with `VolumeIntegralShockCapturingHG` the indicator is
+    # called here to get up-to-date values for IO, this is not easily possible in this case
+    # because the calculation is very integrated into the method.
+    # See also https://github.com/trixi-framework/Trixi.jl/pull/1611#discussion_r1334553206.
+    # Therefore, the coefficients at `t=t^{n-1}` are saved. Thus, the coefficients of the first
+    # stored solution (initial condition) are not yet defined and were manually set to `NaN`.
+    get_node_variables!(node_variables, volume_integral.limiter, volume_integral,
+                        equations)
+end
+
 # TODO: FD. Should this definition live in a different file because it is
 # not strictly a DG method?
 """
@@ -403,6 +420,10 @@ function get_element_variables!(element_variables, u, mesh, equations, dg::DG, c
                            dg, cache)
 end
 
+function get_node_variables!(node_variables, mesh, equations, dg::DG, cache)
+    get_node_variables!(node_variables, mesh, equations, dg.volume_integral, dg, cache)
+end
+
 const MeshesDGSEM = Union{TreeMesh, StructuredMesh, UnstructuredMesh2D, P4estMesh,
                           T8codeMesh}
 

src/solvers/dgsem_tree/containers_2d.jl

@@ -1372,6 +1372,7 @@ function Base.resize!(container::ContainerSubcellLimiterIDP2D, capacity)
     (; _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)
     resize!(_alpha1, (n_nodes + 1) * n_nodes * capacity)
     container.alpha1 = unsafe_wrap(Array, pointer(_alpha1),
                                    (n_nodes + 1, n_nodes, capacity))

src/solvers/dgsem_tree/subcell_limiters.jl

@@ -101,4 +101,11 @@ function Base.show(io::IO, ::MIME"text/plain", limiter::SubcellLimiterIDP)
         summary_box(io, "SubcellLimiterIDP", setup)
     end
 end
+
+function get_node_variables!(node_variables, limiter::SubcellLimiterIDP,
+                             ::VolumeIntegralSubcellLimiting, equations)
+    node_variables[:limiting_coefficient] = limiter.cache.subcell_limiter_coefficients.alpha
+
+    return nothing
+end
 end # @muladd

src/time_integration/methods_SSP.jl

@@ -31,6 +31,11 @@ struct SimpleSSPRK33{StageCallbacks} <: SimpleAlgorithmSSP
     stage_callbacks::StageCallbacks
 
     function SimpleSSPRK33(; stage_callbacks = ())
+        # Mathematically speaking, it is not necessary for the algorithm to split the factors
+        # into numerator and denominator. Otherwise, however, rounding errors of the order of
+        # the machine accuracy will occur, which will add up over time and thus endanger the
+        # conservation of the simulation.
+        # See also https://github.com/trixi-framework/Trixi.jl/pull/1640.
         numerator_a = SVector(0.0, 3.0, 1.0) # a = numerator_a / denominator
         numerator_b = SVector(1.0, 1.0, 2.0) # b = numerator_b / denominator
         denominator = SVector(1.0, 4.0, 3.0)

utils/trixi2txt.jl

@@ -70,7 +70,7 @@ function trixi2txt(filename::AbstractString...;
         center_level_0, length_level_0, leaf_cells, coordinates, levels = read_meshfile(meshfile)
 
         # Read data
-        labels, data, n_elements, n_nodes, element_variables, time = read_datafile(filename)
+        labels, data, n_elements, n_nodes, element_variables, node_variables, time = read_datafile(filename)
 
         # Check if dimensions match
         if length(leaf_cells) != n_elements
@@ -263,7 +263,16 @@ function read_datafile(filename::String)
             index += 1
         end
 
-        return labels, data, n_elements, n_nodes, element_variables, time
+        # Extract node variable arrays
+        node_variables = Dict{String, Union{Vector{Float64}, Vector{Int}}}()
+        index = 1
+        while haskey(file, "node_variables_$index")
+            varname = read(attributes(file["node_variables_$index"])["name"])
+            node_variables[varname] = read(file["node_variables_$index"])
+            index += 1
+        end
+
+        return labels, data, n_elements, n_nodes, element_variables, node_variables, time
     end
 end