Node-level visualization support for coefficients of Subcell limiting

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,12 @@ function Base.show(io::IO, mime::MIME"text/plain",
     end
 end
 
+function get_node_variables!(node_variables, mesh, equations,
+                             volume_integral::VolumeIntegralSubcellLimiting, dg, cache)
+    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 +414,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/subcell_limiters.jl

@@ -100,4 +100,10 @@ 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[:alpha_limiter] = limiter.cache.subcell_limiter_coefficients.alpha
+    return nothing
+end
 end # @muladd

src/solvers/dgsem_tree/subcell_limiters_2d.jl

@@ -6,8 +6,7 @@
 #! format: noindent
 
 # this method is used when the limiter is constructed as for shock-capturing volume integrals
-function create_cache(indicator::Type{SubcellLimiterIDP},
-                      equations::AbstractEquations{2},
+function create_cache(limiter::Type{SubcellLimiterIDP}, equations::AbstractEquations{2},
                       basis::LobattoLegendreBasis, number_bounds)
     subcell_limiter_coefficients = Trixi.ContainerSubcellLimiterIDP2D{real(basis)
                                                                       }(0,

src/time_integration/methods_SSP.jl

@@ -237,5 +237,7 @@ function Base.resize!(semi, volume_integral::VolumeIntegralSubcellLimiting, new_
     # Resize container subcell_limiter_coefficients
     @unpack limiter = volume_integral
     resize!(limiter.cache.subcell_limiter_coefficients, new_size)
+    (; alpha) = limiter.cache.subcell_limiter_coefficients
+    alpha .= zero(eltype(alpha))
 end
 end # @muladd

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