diff --git a/src/callbacks_step/save_solution.jl b/src/callbacks_step/save_solution.jl index 14ea33368f8..31fe0e87c77 100644 --- a/src/callbacks_step/save_solution.jl +++ b/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? diff --git a/src/callbacks_step/save_solution_dg.jl b/src/callbacks_step/save_solution_dg.jl index 6d5004ff65f..7c015999035 100644 --- a/src/callbacks_step/save_solution_dg.jl +++ b/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 diff --git a/src/semidiscretization/semidiscretization.jl b/src/semidiscretization/semidiscretization.jl index c784f716426..fe7858e31ee 100644 --- a/src/semidiscretization/semidiscretization.jl +++ b/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 diff --git a/src/solvers/dg.jl b/src/solvers/dg.jl index 36bbc6de361..30b4f67474f 100644 --- a/src/solvers/dg.jl +++ b/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} diff --git a/src/solvers/dgsem_tree/containers_2d.jl b/src/solvers/dgsem_tree/containers_2d.jl index fc916b41dd2..9e9fe88c15b 100644 --- a/src/solvers/dgsem_tree/containers_2d.jl +++ b/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)) diff --git a/src/solvers/dgsem_tree/subcell_limiters.jl b/src/solvers/dgsem_tree/subcell_limiters.jl index 4d9eec25a89..6197ca6b85d 100644 --- a/src/solvers/dgsem_tree/subcell_limiters.jl +++ b/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 diff --git a/src/time_integration/methods_SSP.jl b/src/time_integration/methods_SSP.jl index 33eb6ebf926..5b72682d48e 100644 --- a/src/time_integration/methods_SSP.jl +++ b/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) diff --git a/utils/trixi2txt.jl b/utils/trixi2txt.jl index b386f150da4..12a3d46760e 100644 --- a/utils/trixi2txt.jl +++ b/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