From 10c8d948a9fce3991bf6f21ac761f513de978a56 Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Tue, 10 Oct 2023 11:59:50 +0200 Subject: [PATCH 1/8] Print addition to DOF count (#1669) --- src/callbacks_step/analysis.jl | 2 +- src/semidiscretization/semidiscretization_coupled.jl | 2 +- src/semidiscretization/semidiscretization_hyperbolic.jl | 2 +- .../semidiscretization_hyperbolic_parabolic.jl | 2 +- 4 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/callbacks_step/analysis.jl b/src/callbacks_step/analysis.jl index fad42b11098..a3e49007fb6 100644 --- a/src/callbacks_step/analysis.jl +++ b/src/callbacks_step/analysis.jl @@ -297,7 +297,7 @@ function (analysis_callback::AnalysisCallback)(u_ode, du_ode, integrator, semi) mpi_println(" " * " " * " " * " PID: " * @sprintf("%10.8e s", performance_index)) - mpi_println(" #DOF: " * @sprintf("% 14d", ndofs(semi)) * + mpi_println(" #DOFs per field:" * @sprintf("% 14d", ndofs(semi)) * " " * " alloc'd memory: " * @sprintf("%14.3f MiB", memory_use)) mpi_println(" #elements: " * diff --git a/src/semidiscretization/semidiscretization_coupled.jl b/src/semidiscretization/semidiscretization_coupled.jl index b7adff78425..49763b12e5d 100644 --- a/src/semidiscretization/semidiscretization_coupled.jl +++ b/src/semidiscretization/semidiscretization_coupled.jl @@ -71,7 +71,7 @@ function Base.show(io::IO, ::MIME"text/plain", semi::SemidiscretizationCoupled) summary_line(increment_indent(io), "source terms", semi.semis[i].source_terms) summary_line(increment_indent(io), "solver", solver |> typeof |> nameof) end - summary_line(io, "total #DOFs", ndofs(semi)) + summary_line(io, "total #DOFs per field", ndofs(semi)) summary_footer(io) end end diff --git a/src/semidiscretization/semidiscretization_hyperbolic.jl b/src/semidiscretization/semidiscretization_hyperbolic.jl index 50b2c21c14e..9d465bfcc5f 100644 --- a/src/semidiscretization/semidiscretization_hyperbolic.jl +++ b/src/semidiscretization/semidiscretization_hyperbolic.jl @@ -243,7 +243,7 @@ function Base.show(io::IO, ::MIME"text/plain", semi::SemidiscretizationHyperboli summary_line(io, "source terms", semi.source_terms) summary_line(io, "solver", semi.solver |> typeof |> nameof) - summary_line(io, "total #DOFs", ndofs(semi)) + summary_line(io, "total #DOFs per field", ndofs(semi)) summary_footer(io) end end diff --git a/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl b/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl index b12ecadb58b..35340d65b1e 100644 --- a/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl +++ b/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl @@ -228,7 +228,7 @@ function Base.show(io::IO, ::MIME"text/plain", summary_line(io, "source terms", semi.source_terms) summary_line(io, "solver", semi.solver |> typeof |> nameof) summary_line(io, "parabolic solver", semi.solver_parabolic |> typeof |> nameof) - summary_line(io, "total #DOFs", ndofs(semi)) + summary_line(io, "total #DOFs per field", ndofs(semi)) summary_footer(io) end end From 8bc5469b3d883be0ce088386e1ec794e448003ef Mon Sep 17 00:00:00 2001 From: Jesse Chan <1156048+jlchan@users.noreply.github.com> Date: Wed, 11 Oct 2023 06:22:49 -0500 Subject: [PATCH 2/8] fix bug in `apply_jacobian_parabolic!` for `P4estMesh` (#1668) * fix bug in parabolic jacobian computations * formatting * updating test results and removing duplicate test --------- Co-authored-by: Hendrik Ranocha --- src/solvers/dgsem_tree/dg_2d_parabolic.jl | 20 +++++++++++++++++++- src/solvers/dgsem_tree/dg_3d_parabolic.jl | 23 +++++++++++++++++++++-- test/test_parabolic_2d.jl | 16 ++++------------ 3 files changed, 44 insertions(+), 15 deletions(-) diff --git a/src/solvers/dgsem_tree/dg_2d_parabolic.jl b/src/solvers/dgsem_tree/dg_2d_parabolic.jl index 3dbc55412ad..1c32703c7c3 100644 --- a/src/solvers/dgsem_tree/dg_2d_parabolic.jl +++ b/src/solvers/dgsem_tree/dg_2d_parabolic.jl @@ -935,7 +935,7 @@ end # This is because the parabolic fluxes are assumed to be of the form # `du/dt + df/dx = dg/dx + source(x,t)`, # where f(u) is the inviscid flux and g(u) is the viscous flux. -function apply_jacobian_parabolic!(du, mesh::Union{TreeMesh{2}, P4estMesh{2}}, +function apply_jacobian_parabolic!(du, mesh::TreeMesh{2}, equations::AbstractEquationsParabolic, dg::DG, cache) @unpack inverse_jacobian = cache.elements @@ -951,4 +951,22 @@ function apply_jacobian_parabolic!(du, mesh::Union{TreeMesh{2}, P4estMesh{2}}, return nothing end + +function apply_jacobian_parabolic!(du, mesh::P4estMesh{2}, + equations::AbstractEquationsParabolic, + dg::DG, cache) + @unpack inverse_jacobian = cache.elements + + @threaded for element in eachelement(dg, cache) + for j in eachnode(dg), i in eachnode(dg) + factor = inverse_jacobian[i, j, element] + + for v in eachvariable(equations) + du[v, i, j, element] *= factor + end + end + end + + return nothing +end end # @muladd diff --git a/src/solvers/dgsem_tree/dg_3d_parabolic.jl b/src/solvers/dgsem_tree/dg_3d_parabolic.jl index 9817e0e5f0e..37492dbcb91 100644 --- a/src/solvers/dgsem_tree/dg_3d_parabolic.jl +++ b/src/solvers/dgsem_tree/dg_3d_parabolic.jl @@ -1125,8 +1125,9 @@ end # This is because the parabolic fluxes are assumed to be of the form # `du/dt + df/dx = dg/dx + source(x,t)`, # where f(u) is the inviscid flux and g(u) is the viscous flux. -function apply_jacobian_parabolic!(du, mesh::Union{TreeMesh{3}, P4estMesh{3}}, - equations::AbstractEquationsParabolic, dg::DG, cache) +function apply_jacobian_parabolic!(du, mesh::TreeMesh{3}, + equations::AbstractEquationsParabolic, + dg::DG, cache) @unpack inverse_jacobian = cache.elements @threaded for element in eachelement(dg, cache) @@ -1141,4 +1142,22 @@ function apply_jacobian_parabolic!(du, mesh::Union{TreeMesh{3}, P4estMesh{3}}, return nothing end + +function apply_jacobian_parabolic!(du, mesh::P4estMesh{3}, + equations::AbstractEquationsParabolic, + dg::DG, cache) + @unpack inverse_jacobian = cache.elements + + @threaded for element in eachelement(dg, cache) + for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) + factor = inverse_jacobian[i, j, k, element] + + for v in eachvariable(equations) + du[v, i, j, k, element] *= factor + end + end + end + + return nothing +end end # @muladd diff --git a/test/test_parabolic_2d.jl b/test/test_parabolic_2d.jl index 967aa1069a0..98aa337afb4 100644 --- a/test/test_parabolic_2d.jl +++ b/test/test_parabolic_2d.jl @@ -260,24 +260,16 @@ isdir(outdir) && rm(outdir, recursive=true) @trixi_testset "P4estMesh2D: elixir_advection_diffusion_periodic_curved.jl" begin @test_trixi_include(joinpath(examples_dir(), "p4est_2d_dgsem", "elixir_advection_diffusion_periodic_curved.jl"), trees_per_dimension = (1, 1), initial_refinement_level = 2, tspan=(0.0, 0.5), - l2 = [0.012380458938507371], - linf = [0.10860506906472567] - ) - end - - @trixi_testset "P4estMesh2D: elixir_advection_diffusion_periodic_curved.jl" begin - @test_trixi_include(joinpath(examples_dir(), "p4est_2d_dgsem", "elixir_advection_diffusion_periodic_curved.jl"), - trees_per_dimension = (1, 1), initial_refinement_level = 2, tspan=(0.0, 0.5), - l2 = [0.012380458938507371], - linf = [0.10860506906472567] + l2 = [0.006708147442490916], + linf = [0.04807038397976693] ) end @trixi_testset "P4estMesh2D: elixir_advection_diffusion_nonperiodic_curved.jl" begin @test_trixi_include(joinpath(examples_dir(), "p4est_2d_dgsem", "elixir_advection_diffusion_nonperiodic_curved.jl"), trees_per_dimension = (1, 1), initial_refinement_level = 2, tspan=(0.0, 0.5), - l2 = [0.04933902988507035], - linf = [0.2550261714590271] + l2 = [0.00919917034843865], + linf = [0.14186297438393505] ) end From 1727b072fa9f269bcde73b516bc7b00f13fc9c99 Mon Sep 17 00:00:00 2001 From: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com> Date: Wed, 11 Oct 2023 16:23:12 +0200 Subject: [PATCH 3/8] Show percentage of simulation progress (#1659) * show percentage of simulation progress * add notes about AnalysisCallback in docs * generalize for general initial time --------- Co-authored-by: Hendrik Ranocha --- docs/src/callbacks.md | 5 +++++ src/callbacks_step/alive.jl | 10 ++++++++-- src/callbacks_step/analysis.jl | 10 ++++++++-- 3 files changed, 21 insertions(+), 4 deletions(-) diff --git a/docs/src/callbacks.md b/docs/src/callbacks.md index 7f44dfd5925..f018bcf7c39 100644 --- a/docs/src/callbacks.md +++ b/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. diff --git a/src/callbacks_step/alive.jl b/src/callbacks_step/alive.jl index eeacd9681d8..9df7181521e 100644 --- a/src/callbacks_step/alive.jl +++ b/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 diff --git a/src/callbacks_step/analysis.jl b/src/callbacks_step/analysis.jl index a3e49007fb6..e5b4a01a885 100644 --- a/src/callbacks_step/analysis.jl +++ b/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(" " * " " * " " * From 7fd4503af9dfdf7022ae75f7398cac8fda811c85 Mon Sep 17 00:00:00 2001 From: Benjamin Bolm <74359358+bennibolm@users.noreply.github.com> Date: Thu, 12 Oct 2023 15:31:43 +0200 Subject: [PATCH 4/8] Subcell limiting: Revise order of bounds using a `Dict` (#1649) * Implement order of bounds as Dict * Implement `variable_bounds` directly as Dictionary * Remove unnecessary line of code * Use Symbols instead of Strings --- src/solvers/dgsem_tree/containers_2d.jl | 31 ++++++++++--------- src/solvers/dgsem_tree/subcell_limiters.jl | 5 +-- src/solvers/dgsem_tree/subcell_limiters_2d.jl | 31 ++++++++----------- 3 files changed, 32 insertions(+), 35 deletions(-) diff --git a/src/solvers/dgsem_tree/containers_2d.jl b/src/solvers/dgsem_tree/containers_2d.jl index 9148b936312..fc916b41dd2 100644 --- a/src/solvers/dgsem_tree/containers_2d.jl +++ b/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)) resize!(_alpha1, (n_nodes + 1) * n_nodes * capacity) @@ -1379,11 +1379,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 diff --git a/src/solvers/dgsem_tree/subcell_limiters.jl b/src/solvers/dgsem_tree/subcell_limiters.jl index 3a707de3bc7..4d9eec25a89 100644 --- a/src/solvers/dgsem_tree/subcell_limiters.jl +++ b/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, diff --git a/src/solvers/dgsem_tree/subcell_limiters_2d.jl b/src/solvers/dgsem_tree/subcell_limiters_2d.jl index 09ab84ed11a..4497217fb56 100644 --- a/src/solvers/dgsem_tree/subcell_limiters_2d.jl +++ b/src/solvers/dgsem_tree/subcell_limiters_2d.jl @@ -6,17 +6,14 @@ #! 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}, - basis::LobattoLegendreBasis, number_bounds) +function create_cache(limiter::Type{SubcellLimiterIDP}, equations::AbstractEquations{2}, + 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, @@ -26,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 @@ -50,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] From dfdaef5adc243fcd6c55262a69ede59fe3a32a8f Mon Sep 17 00:00:00 2001 From: Benjamin Bolm <74359358+bennibolm@users.noreply.github.com> Date: Wed, 18 Oct 2023 16:58:17 +0200 Subject: [PATCH 5/8] 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 | 11 ++++++++-- src/callbacks_step/save_solution_dg.jl | 18 +++++++++++++++-- src/semidiscretization/semidiscretization.jl | 4 ++++ src/solvers/dg.jl | 21 ++++++++++++++++++++ src/solvers/dgsem_tree/containers_2d.jl | 1 + src/solvers/dgsem_tree/subcell_limiters.jl | 7 +++++++ src/time_integration/methods_SSP.jl | 5 +++++ utils/trixi2txt.jl | 13 ++++++++++-- 8 files changed, 74 insertions(+), 6 deletions(-) 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 From bb60b5d973ff4d4e00c37476ef858d37083d2209 Mon Sep 17 00:00:00 2001 From: Simon Candelaresi <10759273+SimonCan@users.noreply.github.com> Date: Thu, 19 Oct 2023 04:51:53 +0100 Subject: [PATCH 6/8] Sc/polytropic (#1526) * Added polytropic Euler equations. * Added polytropic Euler equations. * Added polytropic Euler equations. * Minor clean up in polytropic equations. * Added example for the compressible polytropic Euler equations in 2d on a structured mesh. * Cleaned up example elicir for polytropic Euler in 2d. * Added polytrropic Euler in 2d in tests. * Remoed unused routines. * Update src/equations/polytropic_euler_2d.jl Co-authored-by: Michael Schlottke-Lakemper * Update src/equations/polytropic_euler_2d.jl Co-authored-by: Michael Schlottke-Lakemper * Update src/equations/polytropic_euler_2d.jl Co-authored-by: Michael Schlottke-Lakemper * Update src/equations/polytropic_euler_2d.jl Co-authored-by: Michael Schlottke-Lakemper * Added compressible Euler equations in news. * Auto-reformatted equations. * Auto-reformatted docs for polytropic equations. * Added example elixir for polytropic and isothermal coupling. * Corrected typos. * Update examples/structured_2d_dgsem/elixir_euler_polytropic.jl Co-authored-by: Michael Schlottke-Lakemper * Update NEWS.md Co-authored-by: Michael Schlottke-Lakemper * Renamed elixir_euler_polytropic.jl to elixir_euler_polytropic.jl. * Renamed Euler polytropic elixir in test. * Renamed coupled polytropic example elixir. * Reformatted the boundary conditions for the couple polytropic example elixir. * Removed redundant inv_gamma_minus_one variable. * Update src/equations/polytropic_euler_2d.jl Co-authored-by: Michael Schlottke-Lakemper * Update src/equations/polytropic_euler_2d.jl Co-authored-by: Michael Schlottke-Lakemper * Update src/equations/polytropic_euler_2d.jl Co-authored-by: Michael Schlottke-Lakemper * Update src/equations/polytropic_euler_2d.jl Co-authored-by: Michael Schlottke-Lakemper * Update examples/structured_2d_dgsem/elixir_eulerpolytropic_coupled.jl Yes. That was necessary in an old version for Trixi.jl. Co-authored-by: Michael Schlottke-Lakemper * Update examples/structured_2d_dgsem/elixir_eulerpolytropic_coupled.jl Co-authored-by: Michael Schlottke-Lakemper * Update examples/structured_2d_dgsem/elixir_eulerpolytropic_coupled.jl Co-authored-by: Michael Schlottke-Lakemper * Update examples/structured_2d_dgsem/elixir_eulerpolytropic_coupled.jl Co-authored-by: Michael Schlottke-Lakemper * Removed flux Ranocha for polytropic equations. * Update src/equations/polytropic_euler_2d.jl Co-authored-by: Michael Schlottke-Lakemper * Update src/equations/polytropic_euler_2d.jl Co-authored-by: Hendrik Ranocha * Update src/equations/polytropic_euler_2d.jl Co-authored-by: Hendrik Ranocha * Added Andrew's flux implementation flux_winters_etal. * Added couped polytropic tests. * Added Winter's fluxes for coupled polytropic elixir. * Added flux_winters_etal. * Added convergence test for polytropic Euler. * Renamed polytropic uler convergence Elixir. * Added convergence test initial condition. * Corrected variable names in convergence test for polytropic equations. * Update src/auxiliary/math.jl Co-authored-by: Hendrik Ranocha * Added isotropic Eulre wave test case. * Removed isotropic option in purely polytropic test case. * Implmented Andrew's efficient way of computing coefficients. * Added isotropic Euler wave as test case. * Reformatted polytrropic_2d. * Update examples/structured_2d_dgsem/elixir_eulerpolytropic_wave.jl Co-authored-by: Michael Schlottke-Lakemper * Update src/auxiliary/math.jl Co-authored-by: Andrew Winters * Reformatted example elixir for polytropic/isothermal equations. * Updated errors for elixir_eulerpolytropic_coupled.jl. * Added polytropic convergence test and renamed isothermal file. * Corrected gamma and kappa. * Added source terms for EOC. * Added source term for ECO. * Minor corrections. * Added source terms. * Corrected primary variables. * Removed redundant kappa and gamme redefinitions. * Improved style by using gamma and kappa from equations. * Improved convergence test for polytropic. * Corrected convergence test residual for polytropic equations. * Allow use of flux_lax_friedrichs * Fix comment * Added flux calculation for polytropic Euler. * Update src/equations/polytropic_euler_2d.jl Co-authored-by: Andrew Winters * Update examples/structured_2d_dgsem/elixir_eulerpolytropic_wave.jl Co-authored-by: Andrew Winters * Moved initial condition and source term before flux function. * Update examples/structured_2d_dgsem/elixir_eulerisothermal_wave.jl Co-authored-by: Andrew Winters * Update src/equations/polytropic_euler_2d.jl Co-authored-by: Andrew Winters * Update src/equations/polytropic_euler_2d.jl Co-authored-by: Andrew Winters * Change kappa and gamma in some polytropic Euler elixirs. * Correcte source for Euler polytropic convergence test. * Changed parameters gamma and kappa in Euler polytropic convergence test to other than 1.0. * Reformatted polytropic Euler. * Spelling errors. * Reformatted. * Added elixir for testing entropy conservation for polytropic Euler equations. * Added weak blast wave as initial condition for the polytropic Euler equations. * Added entropy conservation test. * Typo in tests. * Auto reformat. * Corrected polytropic Euler test reults. * Renamed entropy test for Euler polytropic. * Renamed elixir for shock test. * Added coupled Euler polytropic to the tests. * Removed redundant flux function. * Changed example parameters for coupled polytropic Euler. * Removed max-speed calculations, as they are not being used. * Removed coupled polytropic elixir, since this PR is on the polytropic equations only. * Correcte entropy variables for polytropic Euler. This fixes the entropy conservation issue. * Update examples/structured_2d_dgsem/elixir_eulerpolytropic_convergence.jl Co-authored-by: Andrew Winters * Update src/equations/polytropic_euler_2d.jl Co-authored-by: Andrew Winters * Update src/equations/polytropic_euler_2d.jl Co-authored-by: Andrew Winters * Update src/equations/polytropic_euler_2d.jl Co-authored-by: Andrew Winters * Update src/equations/polytropic_euler_2d.jl Co-authored-by: Andrew Winters * Update src/equations/polytropic_euler_2d.jl Co-authored-by: Andrew Winters * Update src/equations/polytropic_euler_2d.jl Co-authored-by: Andrew Winters * Format of comment. * Update src/Trixi.jl Co-authored-by: Michael Schlottke-Lakemper * Updated polytropic wave test results after removing stage limiter. * Removed stqage limiter. * Reformatted Trixi.jl. * Corrected removed exports in Trixi.jl. * Fix for broadcasting. * Renamed elixir for isothermal wave. * Update test/test_structured_2d.jl Co-authored-by: Michael Schlottke-Lakemper * Added comment about isothermal system in elixir. --------- Co-authored-by: Michael Schlottke-Lakemper Co-authored-by: Hendrik Ranocha Co-authored-by: Andrew Winters --- NEWS.md | 1 + .../elixir_eulerpolytropic_convergence.jl | 63 +++++ .../elixir_eulerpolytropic_ec.jl | 80 ++++++ .../elixir_eulerpolytropic_isothermal_wave.jl | 83 ++++++ .../elixir_eulerpolytropic_wave.jl | 80 ++++++ src/Trixi.jl | 5 +- src/auxiliary/math.jl | 67 +++++ src/equations/equations.jl | 5 + src/equations/polytropic_euler_2d.jl | 257 ++++++++++++++++++ test/test_structured_2d.jl | 24 ++ 10 files changed, 663 insertions(+), 2 deletions(-) create mode 100644 examples/structured_2d_dgsem/elixir_eulerpolytropic_convergence.jl create mode 100644 examples/structured_2d_dgsem/elixir_eulerpolytropic_ec.jl create mode 100644 examples/structured_2d_dgsem/elixir_eulerpolytropic_isothermal_wave.jl create mode 100644 examples/structured_2d_dgsem/elixir_eulerpolytropic_wave.jl create mode 100644 src/equations/polytropic_euler_2d.jl diff --git a/NEWS.md b/NEWS.md index 4b96e1e2834..0c78484a782 100644 --- a/NEWS.md +++ b/NEWS.md @@ -12,6 +12,7 @@ for human readability. - Non-uniform `TreeMesh` available for hyperbolic-parabolic equations. - Capability to set truly discontinuous initial conditions in 1D. - Wetting and drying feature and examples for 1D and 2D shallow water equations +- Implementation of the polytropic Euler equations in 2D - Subcell positivity limiting support for conservative variables in 2D for `TreeMesh` #### Changed diff --git a/examples/structured_2d_dgsem/elixir_eulerpolytropic_convergence.jl b/examples/structured_2d_dgsem/elixir_eulerpolytropic_convergence.jl new file mode 100644 index 00000000000..4fc9281e7a0 --- /dev/null +++ b/examples/structured_2d_dgsem/elixir_eulerpolytropic_convergence.jl @@ -0,0 +1,63 @@ + +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the polytropic Euler equations + +gamma = 1.4 +kappa = 0.5 # Scaling factor for the pressure. +equations = PolytropicEulerEquations2D(gamma, kappa) + +initial_condition = initial_condition_convergence_test + +volume_flux = flux_winters_etal +solver = DGSEM(polydeg = 3, surface_flux = flux_hll, + volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) + +coordinates_min = (0.0, 0.0) +coordinates_max = (1.0, 1.0) + +cells_per_dimension = (4, 4) + +mesh = StructuredMesh(cells_per_dimension, + coordinates_min, + coordinates_max) + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + source_terms = source_terms_convergence_test) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 0.1) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 100 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval, + extra_analysis_errors = (:l2_error_primitive, + :linf_error_primitive)) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = 100, + save_initial_solution = true, + save_final_solution = true, + solution_variables = cons2prim) + +stepsize_callback = StepsizeCallback(cfl = 0.1) + +callbacks = CallbackSet(summary_callback, + analysis_callback, alive_callback, + save_solution, + stepsize_callback) + +############################################################################### +# run the simulation + +sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), + dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep = false, callback = callbacks); +summary_callback() # print the timer summary diff --git a/examples/structured_2d_dgsem/elixir_eulerpolytropic_ec.jl b/examples/structured_2d_dgsem/elixir_eulerpolytropic_ec.jl new file mode 100644 index 00000000000..de15f6b2bc3 --- /dev/null +++ b/examples/structured_2d_dgsem/elixir_eulerpolytropic_ec.jl @@ -0,0 +1,80 @@ + +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the polytropic Euler equations + +gamma = 1.4 +kappa = 0.5 # Scaling factor for the pressure. +equations = PolytropicEulerEquations2D(gamma, kappa) + +initial_condition = initial_condition_weak_blast_wave + +############################################################################### +# Get the DG approximation space + +volume_flux = flux_winters_etal +solver = DGSEM(polydeg=4, surface_flux=flux_winters_etal, + volume_integral=VolumeIntegralFluxDifferencing(volume_flux)) + +############################################################################### +# Get the curved quad mesh from a mapping function + +# Mapping as described in https://arxiv.org/abs/2012.12040, but reduced to 2D +function mapping(xi_, eta_) + # Transform input variables between -1 and 1 onto [0,3] + xi = 1.5 * xi_ + 1.5 + eta = 1.5 * eta_ + 1.5 + + y = eta + 3/8 * (cos(1.5 * pi * (2 * xi - 3)/3) * + cos(0.5 * pi * (2 * eta - 3)/3)) + + x = xi + 3/8 * (cos(0.5 * pi * (2 * xi - 3)/3) * + cos(2 * pi * (2 * y - 3)/3)) + + return SVector(x, y) +end + +cells_per_dimension = (16, 16) + +# Create curved mesh with 16 x 16 elements +mesh = StructuredMesh(cells_per_dimension, mapping) + +############################################################################### +# create the semi discretization object + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 2.0) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 100 +analysis_callback = AnalysisCallback(semi, interval=analysis_interval) + +alive_callback = AliveCallback(analysis_interval=analysis_interval) + +save_solution = SaveSolutionCallback(interval=100, + save_initial_solution=true, + save_final_solution=true) + +stepsize_callback = StepsizeCallback(cfl=1.0) + +callbacks = CallbackSet(summary_callback, + analysis_callback, + alive_callback, + save_solution, + stepsize_callback) + +############################################################################### +# run the simulation + +sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false), + dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep=false, callback=callbacks); +summary_callback() # print the timer summary diff --git a/examples/structured_2d_dgsem/elixir_eulerpolytropic_isothermal_wave.jl b/examples/structured_2d_dgsem/elixir_eulerpolytropic_isothermal_wave.jl new file mode 100644 index 00000000000..4ab90957579 --- /dev/null +++ b/examples/structured_2d_dgsem/elixir_eulerpolytropic_isothermal_wave.jl @@ -0,0 +1,83 @@ + +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the polytropic Euler equations + +gamma = 1.0 # With gamma = 1 the system is isothermal. +kappa = 1.0 # Scaling factor for the pressure. +equations = PolytropicEulerEquations2D(gamma, kappa) + +# Linear pressure wave in the negative x-direction. +function initial_condition_wave(x, t, equations::PolytropicEulerEquations2D) + rho = 1.0 + v1 = 0.0 + if x[1] > 0.0 + rho = ((1.0 + 0.01 * sin(x[1] * 2 * pi)) / equations.kappa)^(1 / equations.gamma) + v1 = ((0.01 * sin((x[1] - 1 / 2) * 2 * pi)) / equations.kappa) + end + v2 = 0.0 + + return prim2cons(SVector(rho, v1, v2), equations) +end +initial_condition = initial_condition_wave + +volume_flux = flux_winters_etal +solver = DGSEM(polydeg = 2, surface_flux = flux_hll, + volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) + +coordinates_min = (-2.0, -1.0) +coordinates_max = (2.0, 1.0) + +cells_per_dimension = (64, 32) + +boundary_conditions = (x_neg = boundary_condition_periodic, + x_pos = boundary_condition_periodic, + y_neg = boundary_condition_periodic, + y_pos = boundary_condition_periodic) + +mesh = StructuredMesh(cells_per_dimension, + coordinates_min, + coordinates_max) + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + boundary_conditions = boundary_conditions) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 1.0) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 200 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = 50, + save_initial_solution = true, + save_final_solution = true, + solution_variables = cons2prim) + +stepsize_callback = StepsizeCallback(cfl = 1.7) + +callbacks = CallbackSet(summary_callback, + analysis_callback, alive_callback, + save_solution, + stepsize_callback) + +stage_limiter! = PositivityPreservingLimiterZhangShu(thresholds = (1.0e-4, 1.0e-4), + variables = (Trixi.density, pressure)) + +############################################################################### +# run the simulation + +sol = solve(ode, CarpenterKennedy2N54(stage_limiter!, williamson_condition = false), + dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep = false, callback = callbacks); + +# Print the timer summary +summary_callback() diff --git a/examples/structured_2d_dgsem/elixir_eulerpolytropic_wave.jl b/examples/structured_2d_dgsem/elixir_eulerpolytropic_wave.jl new file mode 100644 index 00000000000..fd332b20aef --- /dev/null +++ b/examples/structured_2d_dgsem/elixir_eulerpolytropic_wave.jl @@ -0,0 +1,80 @@ + +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the polytropic Euler equations + +gamma = 2.0 # Adiabatic monatomic gas in 2d. +kappa = 0.5 # Scaling factor for the pressure. +equations = PolytropicEulerEquations2D(gamma, kappa) + +# Linear pressure wave in the negative x-direction. +function initial_condition_wave(x, t, equations::PolytropicEulerEquations2D) + rho = 1.0 + v1 = 0.0 + if x[1] > 0.0 + rho = ((1.0 + 0.01 * sin(x[1] * 2 * pi)) / equations.kappa)^(1 / equations.gamma) + v1 = ((0.01 * sin((x[1] - 1 / 2) * 2 * pi)) / equations.kappa) + end + v2 = 0.0 + + return prim2cons(SVector(rho, v1, v2), equations) +end +initial_condition = initial_condition_wave + +volume_flux = flux_winters_etal +solver = DGSEM(polydeg = 2, surface_flux = flux_hll, + volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) + +coordinates_min = (-2.0, -1.0) +coordinates_max = (2.0, 1.0) + +cells_per_dimension = (64, 32) + +boundary_conditions = (x_neg = boundary_condition_periodic, + x_pos = boundary_condition_periodic, + y_neg = boundary_condition_periodic, + y_pos = boundary_condition_periodic) + +mesh = StructuredMesh(cells_per_dimension, + coordinates_min, + coordinates_max) + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + boundary_conditions = boundary_conditions) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 1.0) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 200 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = 50, + save_initial_solution = true, + save_final_solution = true, + solution_variables = cons2prim) + +stepsize_callback = StepsizeCallback(cfl = 1.7) + +callbacks = CallbackSet(summary_callback, + analysis_callback, alive_callback, + save_solution, + stepsize_callback) + +############################################################################### +# run the simulation + +sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), + dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep = false, callback = callbacks); + +# Print the timer summary +summary_callback() diff --git a/src/Trixi.jl b/src/Trixi.jl index b65d03e7975..4169700898a 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -150,7 +150,8 @@ export AcousticPerturbationEquations2D, ShallowWaterEquations1D, ShallowWaterEquations2D, ShallowWaterTwoLayerEquations1D, ShallowWaterTwoLayerEquations2D, ShallowWaterEquationsQuasi1D, - LinearizedEulerEquations2D + LinearizedEulerEquations2D, + PolytropicEulerEquations2D export LaplaceDiffusion1D, LaplaceDiffusion2D, CompressibleNavierStokesDiffusion1D, CompressibleNavierStokesDiffusion2D, @@ -165,7 +166,7 @@ export flux, flux_central, flux_lax_friedrichs, flux_hll, flux_hllc, flux_hlle, flux_kennedy_gruber, flux_shima_etal, flux_ec, flux_fjordholm_etal, flux_nonconservative_fjordholm_etal, flux_es_fjordholm_etal, flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal, - flux_chan_etal, flux_nonconservative_chan_etal, + flux_chan_etal, flux_nonconservative_chan_etal, flux_winters_etal, hydrostatic_reconstruction_audusse_etal, flux_nonconservative_audusse_etal, # TODO: TrixiShallowWater: move anything with "chen_noelle" to new file hydrostatic_reconstruction_chen_noelle, flux_nonconservative_chen_noelle, diff --git a/src/auxiliary/math.jl b/src/auxiliary/math.jl index 4ecf7dd3fcc..38ea0bda8c8 100644 --- a/src/auxiliary/math.jl +++ b/src/auxiliary/math.jl @@ -211,4 +211,71 @@ Return `x` if `x` is negative, else zero. In other words, return @inline function negative_part(x) return min(x, zero(x)) end + +""" + stolarsky_mean(x, y, gamma) + +Compute an instance of a weighted Stolarsky mean of the form + + stolarsky_mean(x, y, gamma) = (gamma - 1)/gamma * (y^gamma - x^gamma) / (y^(gamma-1) - x^(gamma-1)) + +where `gamma > 1`. + +Problem: The formula above has a removable singularity at `x == y`. Thus, +some care must be taken to implement it correctly without problems or loss +of accuracy when `x ≈ y`. Here, we use the approach proposed by +Winters et al. (2020). +Set f = (y - x) / (y + x) and g = gamma (for compact notation). +Then, we use the expansions + + ((1+f)^g - (1-f)^g) / g = 2*f + (g-1)(g-2)/3 * f^3 + (g-1)(g-2)(g-3)(g-4)/60 * f^5 + O(f^7) + +and + + ((1+f)^(g-1) - (1-f)^(g-1)) / (g-1) = 2*f + (g-2)(g-3)/3 * f^3 + (g-2)(g-3)(g-4)(g-5)/60 * f^5 + O(f^7) + +Inserting the first few terms of these expansions and performing polynomial long division +we find that + + stolarsky_mean(x, y, gamma) ≈ (y + x) / 2 * (1 + (g-2)/3 * f^2 - (g+1)(g-2)(g-3)/45 * f^4 + (g+1)(g-2)(g-3)(2g(g-2)-9)/945 * f^6) + +Since divisions are usually more expensive on modern hardware than +multiplications (Agner Fog), we try to avoid computing two divisions. Thus, +we use + + f^2 = (y - x)^2 / (x + y)^2 + = (x * (x - 2 * y) + y * y) / (x * (x + 2 * y) + y * y) + +Given ε = 1.0e-4, we use the following algorithm. + + if f^2 < ε + # use the expansion above + else + # use the direct formula (gamma - 1)/gamma * (y^gamma - x^gamma) / (y^(gamma-1) - x^(gamma-1)) + end + +# References +- Andrew R. Winters, Christof Czernik, Moritz B. Schily & Gregor J. Gassner (2020) + Entropy stable numerical approximations for the isothermal and polytropic + Euler equations + [DOI: 10.1007/s10543-019-00789-w](https://doi.org/10.1007/s10543-019-00789-w) +- Agner Fog. + Lists of instruction latencies, throughputs and micro-operation breakdowns + for Intel, AMD, and VIA CPUs. + [https://www.agner.org/optimize/instruction_tables.pdf](https://www.agner.org/optimize/instruction_tables.pdf) +""" +@inline function stolarsky_mean(x, y, gamma) + epsilon_f2 = 1.0e-4 + f2 = (x * (x - 2 * y) + y * y) / (x * (x + 2 * y) + y * y) # f2 = f^2 + if f2 < epsilon_f2 + # convenience coefficients + c1 = (1 / 3) * (gamma - 2) + c2 = -(1 / 15) * (gamma + 1) * (gamma - 3) * c1 + c3 = -(1 / 21) * (2 * gamma * (gamma - 2) - 9) * c2 + return 0.5 * (x + y) * @evalpoly(f2, 1, c1, c2, c3) + else + return (gamma - 1) / gamma * (y^gamma - x^gamma) / + (y^(gamma - 1) - x^(gamma - 1)) + end +end end # @muladd diff --git a/src/equations/equations.jl b/src/equations/equations.jl index 9bae563d85f..3142dcc2765 100644 --- a/src/equations/equations.jl +++ b/src/equations/equations.jl @@ -371,6 +371,11 @@ abstract type AbstractCompressibleEulerMulticomponentEquations{NDIMS, NVARS, NCO include("compressible_euler_multicomponent_1d.jl") include("compressible_euler_multicomponent_2d.jl") +# PolytropicEulerEquations +abstract type AbstractPolytropicEulerEquations{NDIMS, NVARS} <: + AbstractEquations{NDIMS, NVARS} end +include("polytropic_euler_2d.jl") + # Retrieve number of components from equation instance for the multicomponent case @inline function ncomponents(::AbstractCompressibleEulerMulticomponentEquations{NDIMS, NVARS, diff --git a/src/equations/polytropic_euler_2d.jl b/src/equations/polytropic_euler_2d.jl new file mode 100644 index 00000000000..d4902bbafb2 --- /dev/null +++ b/src/equations/polytropic_euler_2d.jl @@ -0,0 +1,257 @@ +# By default, Julia/LLVM does not use fused multiply-add operations (FMAs). +# Since these FMAs can increase the performance of many numerical algorithms, +# we need to opt-in explicitly. +# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details. +@muladd begin +#! format: noindent + +@doc raw""" + PolytropicEulerEquations2D(gamma, kappa) + +The polytropic Euler equations +```math +\frac{\partial}{\partial t} +\begin{pmatrix} +\rho \\ \rho v_1 \\ \rho v_2 +\end{pmatrix} ++ +\frac{\partial}{\partial x} +\begin{pmatrix} + \rho v_1 \\ \rho v_1^2 + \kappa\rho^\gamma \\ \rho v_1 v_2 +\end{pmatrix} ++ +\frac{\partial}{\partial y} +\begin{pmatrix} +\rho v_2 \\ \rho v_1 v_2 \\ \rho v_2^2 + \kappa\rho^\gamma +\end{pmatrix} += +\begin{pmatrix} +0 \\ 0 \\ 0 +\end{pmatrix} +``` +for an ideal gas with ratio of specific heats `gamma` +in two space dimensions. +Here, ``\rho`` is the density and ``v_1`` and`v_2` the velocities and +```math +p = \kappa\rho^\gamma +``` +the pressure, which we replaced using this relation. +""" +struct PolytropicEulerEquations2D{RealT <: Real} <: + AbstractPolytropicEulerEquations{2, 3} + gamma::RealT # ratio of specific heats + kappa::RealT # fluid scaling factor + + function PolytropicEulerEquations2D(gamma, kappa) + gamma_, kappa_ = promote(gamma, kappa) + new{typeof(gamma_)}(gamma_, kappa_) + end +end + +function varnames(::typeof(cons2cons), ::PolytropicEulerEquations2D) + ("rho", "rho_v1", "rho_v2") +end +varnames(::typeof(cons2prim), ::PolytropicEulerEquations2D) = ("rho", "v1", "v2") + +""" + initial_condition_convergence_test(x, t, equations::PolytropicEulerEquations2D) + +Manufactured smooth initial condition used for convergence tests +in combination with [`source_terms_convergence_test`](@ref). +""" +function initial_condition_convergence_test(x, t, equations::PolytropicEulerEquations2D) + # manufactured initial condition from Winters (2019) [0.1007/s10543-019-00789-w] + # domain must be set to [0, 1] x [0, 1] + h = 8 + cos(2 * pi * x[1]) * sin(2 * pi * x[2]) * cos(2 * pi * t) + + return SVector(h, h / 2, 3 * h / 2) +end + +""" + source_terms_convergence_test(u, x, t, equations::PolytropicEulerEquations2D) + +Source terms used for convergence tests in combination with +[`initial_condition_convergence_test`](@ref). +""" +@inline function source_terms_convergence_test(u, x, t, + equations::PolytropicEulerEquations2D) + rho, v1, v2 = cons2prim(u, equations) + + # Residual from Winters (2019) [0.1007/s10543-019-00789-w] eq. (5.2). + h = 8 + cos(2 * pi * x[1]) * sin(2 * pi * x[2]) * cos(2 * pi * t) + h_t = -2 * pi * cos(2 * pi * x[1]) * sin(2 * pi * x[2]) * sin(2 * pi * t) + h_x = -2 * pi * sin(2 * pi * x[1]) * sin(2 * pi * x[2]) * cos(2 * pi * t) + h_y = 2 * pi * cos(2 * pi * x[1]) * cos(2 * pi * x[2]) * cos(2 * pi * t) + + rho_x = h_x + rho_y = h_y + + b = equations.kappa * equations.gamma * h^(equations.gamma - 1) + + r_1 = h_t + h_x / 2 + 3 / 2 * h_y + r_2 = h_t / 2 + h_x / 4 + b * rho_x + 3 / 4 * h_y + r_3 = 3 / 2 * h_t + 3 / 4 * h_x + 9 / 4 * h_y + b * rho_y + + return SVector(r_1, r_2, r_3) +end + +""" + initial_condition_weak_blast_wave(x, t, equations::PolytropicEulerEquations2D) + +A weak blast wave adapted from +- Sebastian Hennemann, Gregor J. Gassner (2020) + A provably entropy stable subcell shock capturing approach for high order split form DG + [arXiv: 2008.12044](https://arxiv.org/abs/2008.12044) +""" +function initial_condition_weak_blast_wave(x, t, equations::PolytropicEulerEquations2D) + # Adapted MHD version of the weak blast wave from Hennemann & Gassner JCP paper 2020 (Sec. 6.3) + # Set up polar coordinates + inicenter = (0, 0) + x_norm = x[1] - inicenter[1] + y_norm = x[2] - inicenter[2] + r = sqrt(x_norm^2 + y_norm^2) + phi = atan(y_norm, x_norm) + + # Calculate primitive variables + rho = r > 0.5 ? 1.0 : 1.1691 + v1 = r > 0.5 ? 0.0 : 0.1882 * cos(phi) + v2 = r > 0.5 ? 0.0 : 0.1882 * sin(phi) + + return prim2cons(SVector(rho, v1, v2), equations) +end + +# Calculate 1D flux for a single point in the normal direction +# Note, this directional vector is not normalized +@inline function flux(u, normal_direction::AbstractVector, + equations::PolytropicEulerEquations2D) + rho, v1, v2 = cons2prim(u, equations) + p = pressure(u, equations) + + v_normal = v1 * normal_direction[1] + v2 * normal_direction[2] + rho_v_normal = rho * v_normal + f1 = rho_v_normal + f2 = rho_v_normal * v1 + p * normal_direction[1] + f3 = rho_v_normal * v2 + p * normal_direction[2] + return SVector(f1, f2, f3) +end + +""" + flux_winters_etal(u_ll, u_rr, normal_direction, + equations::PolytropicEulerEquations2D) + +Entropy conserving two-point flux for isothermal or polytropic gases. +Requires a special weighted Stolarsky mean for the evaluation of the density +denoted here as `stolarsky_mean`. Note, for isothermal gases where `gamma = 1` +this `stolarsky_mean` becomes the [`ln_mean`](@ref). + +For details see Section 3.2 of the following reference +- Andrew R. Winters, Christof Czernik, Moritz B. Schily & Gregor J. Gassner (2020) + Entropy stable numerical approximations for the isothermal and polytropic + Euler equations + [DOI: 10.1007/s10543-019-00789-w](https://doi.org/10.1007/s10543-019-00789-w) +""" +@inline function flux_winters_etal(u_ll, u_rr, normal_direction::AbstractVector, + equations::PolytropicEulerEquations2D) + # Unpack left and right state + rho_ll, v1_ll, v2_ll = cons2prim(u_ll, equations) + rho_rr, v1_rr, v2_rr = cons2prim(u_rr, equations) + p_ll = equations.kappa * rho_ll^equations.gamma + p_rr = equations.kappa * rho_rr^equations.gamma + v_dot_n_ll = v1_ll * normal_direction[1] + v2_ll * normal_direction[2] + v_dot_n_rr = v1_rr * normal_direction[1] + v2_rr * normal_direction[2] + + # Compute the necessary mean values + if equations.gamma == 1.0 # isothermal gas + rho_mean = ln_mean(rho_ll, rho_rr) + else # equations.gamma > 1 # polytropic gas + rho_mean = stolarsky_mean(rho_ll, rho_rr, equations.gamma) + end + v1_avg = 0.5 * (v1_ll + v1_rr) + v2_avg = 0.5 * (v2_ll + v2_rr) + p_avg = 0.5 * (p_ll + p_rr) + + # Calculate fluxes depending on normal_direction + f1 = rho_mean * 0.5 * (v_dot_n_ll + v_dot_n_rr) + f2 = f1 * v1_avg + p_avg * normal_direction[1] + f3 = f1 * v2_avg + p_avg * normal_direction[2] + + return SVector(f1, f2, f3) +end + +@inline function min_max_speed_naive(u_ll, u_rr, normal_direction::AbstractVector, + equations::PolytropicEulerEquations2D) + rho_ll, v1_ll, v2_ll = cons2prim(u_ll, equations) + rho_rr, v1_rr, v2_rr = cons2prim(u_rr, equations) + p_ll = equations.kappa * rho_ll^equations.gamma + p_rr = equations.kappa * rho_rr^equations.gamma + + v_normal_ll = v1_ll * normal_direction[1] + v2_ll * normal_direction[2] + v_normal_rr = v1_rr * normal_direction[1] + v2_rr * normal_direction[2] + + norm_ = norm(normal_direction) + # The v_normals are already scaled by the norm + lambda_min = v_normal_ll - sqrt(equations.gamma * p_ll / rho_ll) * norm_ + lambda_max = v_normal_rr + sqrt(equations.gamma * p_rr / rho_rr) * norm_ + + return lambda_min, lambda_max +end + +@inline function max_abs_speeds(u, equations::PolytropicEulerEquations2D) + rho, v1, v2 = cons2prim(u, equations) + c = sqrt(equations.gamma * equations.kappa * rho^(equations.gamma - 1)) + + return abs(v1) + c, abs(v2) + c +end + +# Convert conservative variables to primitive +@inline function cons2prim(u, equations::PolytropicEulerEquations2D) + rho, rho_v1, rho_v2 = u + + v1 = rho_v1 / rho + v2 = rho_v2 / rho + + return SVector(rho, v1, v2) +end + +# Convert conservative variables to entropy +@inline function cons2entropy(u, equations::PolytropicEulerEquations2D) + rho, rho_v1, rho_v2 = u + + v1 = rho_v1 / rho + v2 = rho_v2 / rho + v_square = v1^2 + v2^2 + p = pressure(u, equations) + # Form of the internal energy depends on gas type + if equations.gamma == 1.0 # isothermal gas + internal_energy = equations.kappa * log(rho) + else # equations.gamma > 1 # polytropic gas + internal_energy = equations.kappa * rho^(equations.gamma - 1) / + (equations.gamma - 1.0) + end + + w1 = internal_energy + p / rho - 0.5 * v_square + w2 = v1 + w3 = v2 + + return SVector(w1, w2, w3) +end + +# Convert primitive to conservative variables +@inline function prim2cons(prim, equations::PolytropicEulerEquations2D) + rho, v1, v2 = prim + rho_v1 = rho * v1 + rho_v2 = rho * v2 + return SVector(rho, rho_v1, rho_v2) +end + +@inline function density(u, equations::PolytropicEulerEquations2D) + rho = u[1] + return rho +end + +@inline function pressure(u, equations::PolytropicEulerEquations2D) + rho, rho_v1, rho_v2 = u + p = equations.kappa * rho^equations.gamma + return p +end +end # @muladd diff --git a/test/test_structured_2d.jl b/test/test_structured_2d.jl index 75937ba82ad..6d528abc7af 100644 --- a/test/test_structured_2d.jl +++ b/test/test_structured_2d.jl @@ -220,6 +220,30 @@ isdir(outdir) && rm(outdir, recursive=true) tspan = (0.0, 0.3)) end + @trixi_testset "elixir_eulerpolytropic_convergence.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulerpolytropic_convergence.jl"), + l2 = [0.0016688820596537988, 0.0025921681885685425, 0.003280950351435014], + linf = [0.010994679664394269, 0.01331197845637, 0.020080117011346488]) + end + + @trixi_testset "elixir_eulerpolytropic_ec.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulerpolytropic_ec.jl"), + l2 = [0.03647890611450939, 0.025284915444045052, 0.025340697771609126], + linf = [0.32516731565355583, 0.37509762516540046, 0.29812843284727336]) + end + + @trixi_testset "elixir_eulerpolytropic_isothermal_wave.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulerpolytropic_isothermal_wave.jl"), + l2 = [0.004998778491726366, 0.004998916000294425, 9.259136963058664e-17], + linf = [0.010001103673834888, 0.010051165098399503, 7.623942913643681e-16]) + end + + @trixi_testset "elixir_eulerpolytropic_wave.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulerpolytropic_wave.jl"), + l2 = [0.23642682112204072, 0.20904264390331334, 8.174982691297391e-17], + linf = [0.4848250368349989, 0.253350873815695, 4.984552457753618e-16]) + end + @trixi_testset "elixir_hypdiff_nonperiodic.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_hypdiff_nonperiodic.jl"), l2 = [0.8799744480157664, 0.8535008397034816, 0.7851383019164209], From fc00ac1bd8a527ac7ed14d9b274b222c874dbf46 Mon Sep 17 00:00:00 2001 From: Benjamin Bolm <74359358+bennibolm@users.noreply.github.com> Date: Fri, 20 Oct 2023 16:16:49 +0200 Subject: [PATCH 7/8] Add implementation of subcell bounds check (#1672) * Add implementation of bounds check * Shorten code by relocate use of `interval` * Add note in docstring * Fix note in docstring of `SubcellLimiterIDP` * Add test for deviations file * Adapt last commit to specific setup in test * Adapt test * Disable test if `!coverage` * Remove non-needed code * Implement first suggestions * Dispatch `init/finalize_callback()` routines * Revise docstring * Rework calculation of bounds * Fix typo * Add comment * Replace construction of `Symbol`s * Implement suggestions for docstring --- .../elixir_euler_shockcapturing_subcell.jl | 2 +- ...ubble_shockcapturing_subcell_positivity.jl | 4 +- src/Trixi.jl | 2 +- src/callbacks_stage/callbacks_stage.jl | 1 + src/callbacks_stage/subcell_bounds_check.jl | 130 ++++++++++++++++++ .../subcell_bounds_check_2d.jl | 49 +++++++ src/solvers/dgsem_tree/subcell_limiters.jl | 4 +- src/solvers/dgsem_tree/subcell_limiters_2d.jl | 12 +- test/test_tree_2d_eulermulti.jl | 7 +- 9 files changed, 203 insertions(+), 8 deletions(-) create mode 100644 src/callbacks_stage/subcell_bounds_check.jl create mode 100644 src/callbacks_stage/subcell_bounds_check_2d.jl diff --git a/examples/tree_2d_dgsem/elixir_euler_shockcapturing_subcell.jl b/examples/tree_2d_dgsem/elixir_euler_shockcapturing_subcell.jl index 6b69e4db563..c696e2de416 100644 --- a/examples/tree_2d_dgsem/elixir_euler_shockcapturing_subcell.jl +++ b/examples/tree_2d_dgsem/elixir_euler_shockcapturing_subcell.jl @@ -84,7 +84,7 @@ callbacks = CallbackSet(summary_callback, ############################################################################### # run the simulation -stage_callbacks = (SubcellLimiterIDPCorrection(),) +stage_callbacks = (SubcellLimiterIDPCorrection(), BoundsCheckCallback(save_errors=false)) sol = Trixi.solve(ode, Trixi.SimpleSSPRK33(stage_callbacks=stage_callbacks); dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback diff --git a/examples/tree_2d_dgsem/elixir_eulermulti_shock_bubble_shockcapturing_subcell_positivity.jl b/examples/tree_2d_dgsem/elixir_eulermulti_shock_bubble_shockcapturing_subcell_positivity.jl index a67eaeb5b2b..c5a7a5932e6 100644 --- a/examples/tree_2d_dgsem/elixir_eulermulti_shock_bubble_shockcapturing_subcell_positivity.jl +++ b/examples/tree_2d_dgsem/elixir_eulermulti_shock_bubble_shockcapturing_subcell_positivity.jl @@ -132,7 +132,9 @@ callbacks = CallbackSet(summary_callback, ############################################################################### # run the simulation -stage_callbacks = (SubcellLimiterIDPCorrection(),) +output_directory = "out" +stage_callbacks = (SubcellLimiterIDPCorrection(), + BoundsCheckCallback(save_errors=true, interval=100, output_directory=output_directory)) sol = Trixi.solve(ode, Trixi.SimpleSSPRK33(stage_callbacks=stage_callbacks); dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback diff --git a/src/Trixi.jl b/src/Trixi.jl index 4169700898a..457d9dc336d 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -232,7 +232,7 @@ export DG, SurfaceIntegralUpwind, MortarL2 -export VolumeIntegralSubcellLimiting, +export VolumeIntegralSubcellLimiting, BoundsCheckCallback, SubcellLimiterIDP, SubcellLimiterIDPCorrection export nelements, nnodes, nvariables, diff --git a/src/callbacks_stage/callbacks_stage.jl b/src/callbacks_stage/callbacks_stage.jl index 976af327e6f..70d60de7914 100644 --- a/src/callbacks_stage/callbacks_stage.jl +++ b/src/callbacks_stage/callbacks_stage.jl @@ -7,6 +7,7 @@ include("positivity_zhang_shu.jl") include("subcell_limiter_idp_correction.jl") +include("subcell_bounds_check.jl") # TODO: TrixiShallowWater: move specific limiter file include("positivity_shallow_water.jl") end # @muladd diff --git a/src/callbacks_stage/subcell_bounds_check.jl b/src/callbacks_stage/subcell_bounds_check.jl new file mode 100644 index 00000000000..c86f266147c --- /dev/null +++ b/src/callbacks_stage/subcell_bounds_check.jl @@ -0,0 +1,130 @@ +# By default, Julia/LLVM does not use fused multiply-add operations (FMAs). +# Since these FMAs can increase the performance of many numerical algorithms, +# we need to opt-in explicitly. +# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details. +@muladd begin +#! format: noindent + +""" + BoundsCheckCallback(; output_directory="out", save_errors=false, interval=1) + +Subcell limiting techniques with [`SubcellLimiterIDP`](@ref) are constructed to adhere certain +local or global bounds. To make sure that these bounds are actually met, this callback calculates +the maximum deviation from the bounds. The maximum deviation per applied bound is printed to +the screen at the end of the simulation. +For more insights, when setting `save_errors=true` the occurring errors are exported every +`interval` time steps during the simulation. Then, the maximum deviations since the last +export are saved in "`output_directory`/deviations.txt". +The `BoundsCheckCallback` has to be applied as a stage callback for the SSPRK time integration scheme. + +!!! note + For `SubcellLimiterIDP`, the solution is corrected in the a posteriori correction stage + [`SubcellLimiterIDPCorrection`](@ref). So, to check the final solution, this bounds check + callback must be called after the correction stage. +""" +struct BoundsCheckCallback + output_directory::String + save_errors::Bool + interval::Int +end + +function BoundsCheckCallback(; output_directory = "out", save_errors = false, + interval = 1) + BoundsCheckCallback(output_directory, save_errors, interval) +end + +function (callback::BoundsCheckCallback)(u_ode, integrator, stage) + mesh, equations, solver, cache = mesh_equations_solver_cache(integrator.p) + (; t, iter, alg) = integrator + u = wrap_array(u_ode, mesh, equations, solver, cache) + + save_errors = callback.save_errors && (callback.interval > 0) && + (stage == length(alg.c)) && + (iter % callback.interval == 0 || integrator.finalstep) + @trixi_timeit timer() "check_bounds" check_bounds(u, mesh, equations, solver, cache, + solver.volume_integral, t, + iter + 1, + callback.output_directory, + save_errors) +end + +function check_bounds(u, mesh, equations, solver, cache, + volume_integral::AbstractVolumeIntegral, t, iter, + output_directory, save_errors) + return nothing +end + +function check_bounds(u, mesh, equations, solver, cache, + volume_integral::VolumeIntegralSubcellLimiting, t, iter, + output_directory, save_errors) + check_bounds(u, mesh, equations, solver, cache, volume_integral.limiter, t, iter, + output_directory, save_errors) +end + +function init_callback(callback::BoundsCheckCallback, semi) + init_callback(callback, semi, semi.solver.volume_integral) +end + +init_callback(callback::BoundsCheckCallback, semi, volume_integral::AbstractVolumeIntegral) = nothing + +function init_callback(callback::BoundsCheckCallback, semi, + volume_integral::VolumeIntegralSubcellLimiting) + init_callback(callback, semi, volume_integral.limiter) +end + +function init_callback(callback::BoundsCheckCallback, semi, limiter::SubcellLimiterIDP) + if !callback.save_errors || (callback.interval == 0) + return nothing + end + + (; positivity) = limiter + (; output_directory) = callback + variables = varnames(cons2cons, semi.equations) + + mkpath(output_directory) + open("$output_directory/deviations.txt", "a") do f + print(f, "# iter, simu_time") + if positivity + for v in limiter.positivity_variables_cons + print(f, ", " * string(variables[v]) * "_min") + end + end + println(f) + end + + return nothing +end + +function finalize_callback(callback::BoundsCheckCallback, semi) + finalize_callback(callback, semi, semi.solver.volume_integral) +end + +finalize_callback(callback::BoundsCheckCallback, semi, volume_integral::AbstractVolumeIntegral) = nothing + +function finalize_callback(callback::BoundsCheckCallback, semi, + volume_integral::VolumeIntegralSubcellLimiting) + finalize_callback(callback, semi, volume_integral.limiter) +end + +@inline function finalize_callback(callback::BoundsCheckCallback, semi, + limiter::SubcellLimiterIDP) + (; positivity) = limiter + (; idp_bounds_delta) = limiter.cache + variables = varnames(cons2cons, semi.equations) + + println("─"^100) + println("Maximum deviation from bounds:") + println("─"^100) + if positivity + for v in limiter.positivity_variables_cons + println(string(variables[v]) * ":\n- positivity: ", + idp_bounds_delta[Symbol(string(v), "_min")][2]) + end + end + println("─"^100 * "\n") + + return nothing +end + +include("subcell_bounds_check_2d.jl") +end # @muladd diff --git a/src/callbacks_stage/subcell_bounds_check_2d.jl b/src/callbacks_stage/subcell_bounds_check_2d.jl new file mode 100644 index 00000000000..8159becb503 --- /dev/null +++ b/src/callbacks_stage/subcell_bounds_check_2d.jl @@ -0,0 +1,49 @@ +# By default, Julia/LLVM does not use fused multiply-add operations (FMAs). +# Since these FMAs can increase the performance of many numerical algorithms, +# we need to opt-in explicitly. +# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details. +@muladd begin +#! format: noindent + +@inline function check_bounds(u, mesh::AbstractMesh{2}, equations, solver, cache, + limiter::SubcellLimiterIDP, + time, iter, output_directory, save_errors) + (; positivity) = solver.volume_integral.limiter + (; variable_bounds) = limiter.cache.subcell_limiter_coefficients + (; idp_bounds_delta) = limiter.cache + + if positivity + for v in limiter.positivity_variables_cons + key = Symbol(string(v), "_min") + deviation = idp_bounds_delta[key] + for element in eachelement(solver, cache), j in eachnode(solver), + i in eachnode(solver) + + var = u[v, i, j, element] + deviation[1] = max(deviation[1], + variable_bounds[key][i, j, element] - var) + end + deviation[2] = max(deviation[2], deviation[1]) + end + end + if save_errors + # Print to output file + open("$output_directory/deviations.txt", "a") do f + print(f, iter, ", ", time) + if positivity + for v in limiter.positivity_variables_cons + key = Symbol(string(v), "_min") + print(f, ", ", idp_bounds_delta[key][1]) + end + end + println(f) + end + # Reset first entries of idp_bounds_delta + for (key, _) in idp_bounds_delta + idp_bounds_delta[key][1] = zero(eltype(idp_bounds_delta[key][1])) + end + end + + return nothing +end +end # @muladd diff --git a/src/solvers/dgsem_tree/subcell_limiters.jl b/src/solvers/dgsem_tree/subcell_limiters.jl index 6197ca6b85d..55d402954bf 100644 --- a/src/solvers/dgsem_tree/subcell_limiters.jl +++ b/src/solvers/dgsem_tree/subcell_limiters.jl @@ -26,7 +26,7 @@ The bounds are calculated using the low-order FV solution. The positivity limite !!! note This limiter and the correction callback [`SubcellLimiterIDPCorrection`](@ref) only work together. - Without the callback, no limiting takes place, leading to a standard flux-differencing DGSEM scheme. + Without the callback, no correction takes place, leading to a standard low-order FV scheme. ## References @@ -53,7 +53,7 @@ function SubcellLimiterIDP(equations::AbstractEquations, basis; positivity_correction_factor = 0.1) positivity = (length(positivity_variables_cons) > 0) - bound_keys = Tuple(Symbol("$(i)_min") for i in positivity_variables_cons) + bound_keys = Tuple(Symbol(string(v), "_min") for v in positivity_variables_cons) cache = create_cache(SubcellLimiterIDP, equations, basis, bound_keys) diff --git a/src/solvers/dgsem_tree/subcell_limiters_2d.jl b/src/solvers/dgsem_tree/subcell_limiters_2d.jl index 4497217fb56..5e00ab4e903 100644 --- a/src/solvers/dgsem_tree/subcell_limiters_2d.jl +++ b/src/solvers/dgsem_tree/subcell_limiters_2d.jl @@ -13,7 +13,15 @@ function create_cache(limiter::Type{SubcellLimiterIDP}, equations::AbstractEquat nnodes(basis), bound_keys) - return (; subcell_limiter_coefficients) + # Memory for bounds checking routine with `BoundsCheckCallback`. + # The first entry of each vector contains the maximum deviation since the last export. + # The second one contains the total maximum deviation. + idp_bounds_delta = Dict{Symbol, Vector{real(basis)}}() + for key in bound_keys + idp_bounds_delta[key] = zeros(real(basis), 2) + end + + return (; subcell_limiter_coefficients, idp_bounds_delta) end function (limiter::SubcellLimiterIDP)(u::AbstractArray{<:Any, 4}, semi, dg::DGSEM, t, @@ -60,7 +68,7 @@ end (; positivity_correction_factor) = limiter (; variable_bounds) = limiter.cache.subcell_limiter_coefficients - var_min = variable_bounds[Symbol("$(variable)_min")] + var_min = variable_bounds[Symbol(string(variable), "_min")] @threaded for element in eachelement(dg, semi.cache) inverse_jacobian = cache.elements.inverse_jacobian[element] diff --git a/test/test_tree_2d_eulermulti.jl b/test/test_tree_2d_eulermulti.jl index 606afca1034..0ae46a92ef8 100644 --- a/test/test_tree_2d_eulermulti.jl +++ b/test/test_tree_2d_eulermulti.jl @@ -20,11 +20,16 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_2d_dgsem") end @trixi_testset "elixir_eulermulti_shock_bubble_shockcapturing_subcell_positivity.jl" begin + rm("out/deviations.txt", force=true) @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulermulti_shock_bubble_shockcapturing_subcell_positivity.jl"), l2 = [81.52845664909304, 2.5455678559421346, 63229.190712645846, 0.19929478404550321, 0.011068604228443425], linf = [249.21708417382013, 40.33299887640794, 174205.0118831558, 0.6881458768113586, 0.11274401158173972], initial_refinement_level = 3, - tspan = (0.0, 0.001)) + tspan = (0.0, 0.001), + output_directory="out") + lines = readlines("out/deviations.txt") + @test lines[1] == "# iter, simu_time, rho1_min, rho2_min" + @test startswith(lines[end], "1") end @trixi_testset "elixir_eulermulti_ec.jl" begin From fcab6a9e50531e5a77e0de77f2f369d75c5a7716 Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Sat, 21 Oct 2023 17:58:50 +0200 Subject: [PATCH 8/8] Add comments explaining why weak form is only for conserved terms (#1676) * Add comments explaining why weak form is only for conserved terms * Update src/solvers/dg.jl Co-authored-by: Andrew Winters * Update src/solvers/dgsem_structured/dg_2d.jl Co-authored-by: Andrew Winters * Update src/solvers/dg.jl Co-authored-by: Andrew Winters * Update src/solvers/dgsem_structured/dg_3d.jl Co-authored-by: Andrew Winters * Update src/solvers/dgsem_structured/dg_2d.jl Co-authored-by: Andrew Winters * Update src/solvers/dgsem_tree/dg_2d.jl Co-authored-by: Andrew Winters * Update src/solvers/dgsem_tree/dg_1d.jl Co-authored-by: Andrew Winters * Update src/solvers/dgsem_tree/dg_3d.jl Co-authored-by: Andrew Winters * Update src/solvers/dgsem_tree/dg_2d.jl Co-authored-by: Andrew Winters * Update src/solvers/dgsem_tree/dg_3d.jl Co-authored-by: Andrew Winters * Update src/solvers/dgsem_structured/dg_3d.jl Co-authored-by: Andrew Winters * Update src/solvers/dgsem_tree/dg_1d.jl Co-authored-by: Andrew Winters * Update src/solvers/dg.jl Co-authored-by: Hendrik Ranocha --------- Co-authored-by: Andrew Winters Co-authored-by: Hendrik Ranocha --- src/solvers/dg.jl | 5 +++++ src/solvers/dgsem_structured/dg_2d.jl | 7 +++++++ src/solvers/dgsem_structured/dg_3d.jl | 7 +++++++ src/solvers/dgsem_tree/dg_1d.jl | 7 +++++++ src/solvers/dgsem_tree/dg_2d.jl | 7 +++++++ src/solvers/dgsem_tree/dg_3d.jl | 7 +++++++ 6 files changed, 40 insertions(+) diff --git a/src/solvers/dg.jl b/src/solvers/dg.jl index 30b4f67474f..91ad59b76b6 100644 --- a/src/solvers/dg.jl +++ b/src/solvers/dg.jl @@ -40,6 +40,11 @@ standard textbooks. Nodal Discontinuous Galerkin Methods: Algorithms, Analysis, and Applications [doi: 10.1007/978-0-387-72067-8](https://doi.org/10.1007/978-0-387-72067-8) + +`VolumeIntegralWeakForm()` is only implemented for conserved terms as +non-conservative terms should always be discretized in conjunction with a flux-splitting scheme, +see [`VolumeIntegralFluxDifferencing`](@ref). +This treatment is required to achieve, e.g., entropy-stability or well-balancedness. """ struct VolumeIntegralWeakForm <: AbstractVolumeIntegral end diff --git a/src/solvers/dgsem_structured/dg_2d.jl b/src/solvers/dgsem_structured/dg_2d.jl index 3e8ce759b30..25a0eea096f 100644 --- a/src/solvers/dgsem_structured/dg_2d.jl +++ b/src/solvers/dgsem_structured/dg_2d.jl @@ -49,6 +49,13 @@ function rhs!(du, u, t, return nothing end +#= +`weak_form_kernel!` is only implemented for conserved terms as +non-conservative terms should always be discretized in conjunction with a flux-splitting scheme, +see `flux_differencing_kernel!`. +This treatment is required to achieve, e.g., entropy-stability or well-balancedness. +See also https://github.com/trixi-framework/Trixi.jl/issues/1671#issuecomment-1765644064 +=# @inline function weak_form_kernel!(du, u, element, mesh::Union{StructuredMesh{2}, UnstructuredMesh2D, diff --git a/src/solvers/dgsem_structured/dg_3d.jl b/src/solvers/dgsem_structured/dg_3d.jl index 0e6bf8a2ac0..cdb085e9008 100644 --- a/src/solvers/dgsem_structured/dg_3d.jl +++ b/src/solvers/dgsem_structured/dg_3d.jl @@ -49,6 +49,13 @@ function rhs!(du, u, t, return nothing end +#= +`weak_form_kernel!` is only implemented for conserved terms as +non-conservative terms should always be discretized in conjunction with a flux-splitting scheme, +see `flux_differencing_kernel!`. +This treatment is required to achieve, e.g., entropy-stability or well-balancedness. +See also https://github.com/trixi-framework/Trixi.jl/issues/1671#issuecomment-1765644064 +=# @inline function weak_form_kernel!(du, u, element, mesh::Union{StructuredMesh{3}, P4estMesh{3}}, diff --git a/src/solvers/dgsem_tree/dg_1d.jl b/src/solvers/dgsem_tree/dg_1d.jl index b5bb076f3b7..4a0747d1c09 100644 --- a/src/solvers/dgsem_tree/dg_1d.jl +++ b/src/solvers/dgsem_tree/dg_1d.jl @@ -144,6 +144,13 @@ function calc_volume_integral!(du, u, return nothing end +#= +`weak_form_kernel!` is only implemented for conserved terms as +non-conservative terms should always be discretized in conjunction with a flux-splitting scheme, +see `flux_differencing_kernel!`. +This treatment is required to achieve, e.g., entropy-stability or well-balancedness. +See also https://github.com/trixi-framework/Trixi.jl/issues/1671#issuecomment-1765644064 +=# @inline function weak_form_kernel!(du, u, element, mesh::Union{TreeMesh{1}, StructuredMesh{1}}, nonconservative_terms::False, equations, diff --git a/src/solvers/dgsem_tree/dg_2d.jl b/src/solvers/dgsem_tree/dg_2d.jl index c30d0a8e01a..7ecf4c00032 100644 --- a/src/solvers/dgsem_tree/dg_2d.jl +++ b/src/solvers/dgsem_tree/dg_2d.jl @@ -194,6 +194,13 @@ function calc_volume_integral!(du, u, return nothing end +#= +`weak_form_kernel!` is only implemented for conserved terms as +non-conservative terms should always be discretized in conjunction with a flux-splitting scheme, +see `flux_differencing_kernel!`. +This treatment is required to achieve, e.g., entropy-stability or well-balancedness. +See also https://github.com/trixi-framework/Trixi.jl/issues/1671#issuecomment-1765644064 +=# @inline function weak_form_kernel!(du, u, element, mesh::TreeMesh{2}, nonconservative_terms::False, equations, diff --git a/src/solvers/dgsem_tree/dg_3d.jl b/src/solvers/dgsem_tree/dg_3d.jl index acdab900cd1..3364187e93c 100644 --- a/src/solvers/dgsem_tree/dg_3d.jl +++ b/src/solvers/dgsem_tree/dg_3d.jl @@ -223,6 +223,13 @@ function calc_volume_integral!(du, u, return nothing end +#= +`weak_form_kernel!` is only implemented for conserved terms as +non-conservative terms should always be discretized in conjunction with a flux-splitting scheme, +see `flux_differencing_kernel!`. +This treatment is required to achieve, e.g., entropy-stability or well-balancedness. +See also https://github.com/trixi-framework/Trixi.jl/issues/1671#issuecomment-1765644064 +=# @inline function weak_form_kernel!(du, u, element, mesh::TreeMesh{3}, nonconservative_terms::False, equations,