From 062568c2dfa826235e3ef1e39c178388eb2cd0f1 Mon Sep 17 00:00:00 2001 From: bennibolm Date: Tue, 14 Nov 2023 16:24:55 +0100 Subject: [PATCH 01/25] Add positivity limiting of non-linear variables --- ...kelvin_helmholtz_instability_sc_subcell.jl | 91 +++++++++ .../elixir_mhd_shockcapturing_subcell.jl | 7 +- src/callbacks_stage/subcell_bounds_check.jl | 8 + .../subcell_bounds_check_2d.jl | 17 ++ src/equations/compressible_euler_2d.jl | 19 ++ src/equations/ideal_glm_mhd_2d.jl | 21 ++ src/solvers/dgsem_tree/subcell_limiters.jl | 43 ++++- src/solvers/dgsem_tree/subcell_limiters_2d.jl | 181 ++++++++++++++++++ test/test_tree_2d_euler.jl | 26 +++ test/test_tree_2d_mhd.jl | 32 ++-- test/test_unit.jl | 3 +- 11 files changed, 422 insertions(+), 26 deletions(-) create mode 100644 examples/tree_2d_dgsem/elixir_euler_kelvin_helmholtz_instability_sc_subcell.jl diff --git a/examples/tree_2d_dgsem/elixir_euler_kelvin_helmholtz_instability_sc_subcell.jl b/examples/tree_2d_dgsem/elixir_euler_kelvin_helmholtz_instability_sc_subcell.jl new file mode 100644 index 00000000000..1817672778a --- /dev/null +++ b/examples/tree_2d_dgsem/elixir_euler_kelvin_helmholtz_instability_sc_subcell.jl @@ -0,0 +1,91 @@ + +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the compressible Euler equations +gamma = 1.4 +equations = CompressibleEulerEquations2D(gamma) + +""" + initial_condition_kelvin_helmholtz_instability(x, t, equations::CompressibleEulerEquations2D) + +A version of the classical Kelvin-Helmholtz instability based on +- Andrés M. Rueda-Ramírez, Gregor J. Gassner (2021) + A Subcell Finite Volume Positivity-Preserving Limiter for DGSEM Discretizations + of the Euler Equations + [arXiv: 2102.06017](https://arxiv.org/abs/2102.06017) +""" +function initial_condition_kelvin_helmholtz_instability(x, t, + equations::CompressibleEulerEquations2D) + # change discontinuity to tanh + # typical resolution 128^2, 256^2 + # domain size is [-1,+1]^2 + slope = 15 + amplitude = 0.02 + B = tanh(slope * x[2] + 7.5) - tanh(slope * x[2] - 7.5) + rho = 0.5 + 0.75 * B + v1 = 0.5 * (B - 1) + v2 = 0.1 * sin(2 * pi * x[1]) + p = 1.0 + return prim2cons(SVector(rho, v1, v2, p), equations) +end +initial_condition = initial_condition_kelvin_helmholtz_instability + +surface_flux = flux_lax_friedrichs +volume_flux = flux_ranocha +polydeg = 3 +basis = LobattoLegendreBasis(polydeg) + +limiter_idp = SubcellLimiterIDP(equations, basis; + positivity_variables_cons = ["rho"], + positivity_variables_nonlinear = [pressure]) +volume_integral = VolumeIntegralSubcellLimiting(limiter_idp; + volume_flux_dg = volume_flux, + volume_flux_fv = surface_flux) +solver = DGSEM(basis, surface_flux, volume_integral) + +coordinates_min = (-1.0, -1.0) +coordinates_max = (1.0, 1.0) +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level = 5, + n_cells_max = 100_000) +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 3.7) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 1000 +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, + solution_variables = cons2prim) + +save_restart = SaveRestartCallback(interval = 1000, + save_final_restart = true) + +stepsize_callback = StepsizeCallback(cfl = 0.7) + +callbacks = CallbackSet(summary_callback, + analysis_callback, alive_callback, + stepsize_callback, + save_restart, save_solution) + +############################################################################### +# run the simulation + +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 + callback = callbacks); +summary_callback() # print the timer summary diff --git a/examples/tree_2d_dgsem/elixir_mhd_shockcapturing_subcell.jl b/examples/tree_2d_dgsem/elixir_mhd_shockcapturing_subcell.jl index fe9ad92467f..74d0370647a 100644 --- a/examples/tree_2d_dgsem/elixir_mhd_shockcapturing_subcell.jl +++ b/examples/tree_2d_dgsem/elixir_mhd_shockcapturing_subcell.jl @@ -22,7 +22,7 @@ function initial_condition_blast_wave(x, t, equations::IdealGlmMhdEquations2D) r = sqrt(x[1]^2 + x[2]^2) pmax = 10.0 - pmin = 1.0 + pmin = 0.01 rhomax = 1.0 rhomin = 0.01 if r <= 0.09 @@ -52,7 +52,8 @@ basis = LobattoLegendreBasis(3) limiter_idp = SubcellLimiterIDP(equations, basis; positivity_variables_cons = ["rho"], - positivity_correction_factor = 0.5) + positivity_variables_nonlinear = [pressure], + positivity_correction_factor = 0.1) volume_integral = VolumeIntegralSubcellLimiting(limiter_idp; volume_flux_dg = volume_flux, volume_flux_fv = surface_flux) @@ -84,7 +85,7 @@ save_solution = SaveSolutionCallback(interval = 100, save_final_solution = true, solution_variables = cons2prim) -cfl = 0.5 +cfl = 0.4 stepsize_callback = StepsizeCallback(cfl = cfl) glm_speed_callback = GlmSpeedCallback(glm_scale = 0.5, cfl = cfl) diff --git a/src/callbacks_stage/subcell_bounds_check.jl b/src/callbacks_stage/subcell_bounds_check.jl index d7e30ab1621..25662f85802 100644 --- a/src/callbacks_stage/subcell_bounds_check.jl +++ b/src/callbacks_stage/subcell_bounds_check.jl @@ -97,6 +97,9 @@ function init_callback(callback::BoundsCheckCallback, semi, limiter::SubcellLimi end print(f, ", " * string(variables[v]) * "_min") end + for variable in limiter.positivity_variables_nonlinear + print(f, ", " * string(variable) * "_min") + end end println(f) end @@ -140,6 +143,11 @@ end println(string(variables[v]) * ":\n- positivity: ", idp_bounds_delta[Symbol(string(v), "_min")][2]) end + for variable in limiter.positivity_variables_nonlinear + variable_string = string(variable) + println(variable_string * ":\n- positivity: ", + idp_bounds_delta[Symbol(variable_string, "_min")][2]) + end end println("─"^100 * "\n") diff --git a/src/callbacks_stage/subcell_bounds_check_2d.jl b/src/callbacks_stage/subcell_bounds_check_2d.jl index d52eb6edb9e..688d3281b9e 100644 --- a/src/callbacks_stage/subcell_bounds_check_2d.jl +++ b/src/callbacks_stage/subcell_bounds_check_2d.jl @@ -48,6 +48,19 @@ end deviation[2] = max(deviation[2], deviation[1]) end + for variable in limiter.positivity_variables_nonlinear + key = Symbol(string(variable), "_min") + deviation = idp_bounds_delta[key] + for element in eachelement(solver, cache), j in eachnode(solver), + i in eachnode(solver) + + var = variable(get_node_vars(u, equations, solver, i, j, element), + equations) + 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 @@ -67,6 +80,10 @@ end print(f, ", ", idp_bounds_delta[Symbol(string(v), "_min")][1]) end + for variable in limiter.positivity_variables_nonlinear + print(f, ", ", + idp_bounds_delta[Symbol(string(variable), "_min")][1]) + end end println(f) end diff --git a/src/equations/compressible_euler_2d.jl b/src/equations/compressible_euler_2d.jl index a992f99eaf4..f3d46da2451 100644 --- a/src/equations/compressible_euler_2d.jl +++ b/src/equations/compressible_euler_2d.jl @@ -1404,6 +1404,17 @@ end return SVector(w1, w2, w3, w4) end +# Transformation from conservative variables u to d(p)/d(u) +@inline function pressure(u, equations::CompressibleEulerEquations2D, derivative::True) + rho, rho_v1, rho_v2, rho_e = u + + v1 = rho_v1 / rho + v2 = rho_v2 / rho + v_square = v1^2 + v2^2 + + return (equations.gamma - 1.0) * SVector(0.5 * v_square, -v1, -v2, 1.0) +end + @inline function entropy2cons(w, equations::CompressibleEulerEquations2D) # See Hughes, Franca, Mallet (1986) A new finite element formulation for CFD # [DOI: 10.1016/0045-7825(86)90127-1](https://doi.org/10.1016/0045-7825(86)90127-1) @@ -1428,6 +1439,14 @@ end return SVector(rho, rho_v1, rho_v2, rho_e) end +@inline function is_valid_state(cons, equations::CompressibleEulerEquations2D) + p = pressure(cons, equations) + if cons[1] <= 0.0 || p <= 0.0 + return false + end + return true +end + # Convert primitive to conservative variables @inline function prim2cons(prim, equations::CompressibleEulerEquations2D) rho, v1, v2, p = prim diff --git a/src/equations/ideal_glm_mhd_2d.jl b/src/equations/ideal_glm_mhd_2d.jl index 43d1991e34b..dd80803e6a6 100644 --- a/src/equations/ideal_glm_mhd_2d.jl +++ b/src/equations/ideal_glm_mhd_2d.jl @@ -1062,6 +1062,19 @@ end return SVector(w1, w2, w3, w4, w5, w6, w7, w8, w9) end +# Transformation from conservative variables u to d(p)/d(u) +@inline function pressure(u, equations::IdealGlmMhdEquations2D, derivative::True) + rho, rho_v1, rho_v2, rho_v3, rho_e, B1, B2, B3, psi = u + + v1 = rho_v1 / rho + v2 = rho_v2 / rho + v3 = rho_v3 / rho + v_square = v1^2 + v2^2 + v3^2 + + return (equations.gamma - 1.0) * + SVector(0.5 * v_square, -v1, -v2, -v3, 1.0, -B1, -B2, -B3, -psi) +end + # Convert entropy variables to conservative variables @inline function entropy2cons(w, equations::IdealGlmMhdEquations2D) w1, w2, w3, w4, w5, w6, w7, w8, w9 = w @@ -1089,6 +1102,14 @@ end return prim2cons(SVector(rho, v1, v2, v3, p, B1, B2, B3, psi), equations) end +@inline function is_valid_state(cons, equations::IdealGlmMhdEquations2D) + p = pressure(cons, equations) + if cons[1] <= 0.0 || p <= 0.0 + return false + end + return true +end + # Convert primitive to conservative variables @inline function prim2cons(prim, equations::IdealGlmMhdEquations2D) rho, v1, v2, v3, p, B1, B2, B3, psi = prim diff --git a/src/solvers/dgsem_tree/subcell_limiters.jl b/src/solvers/dgsem_tree/subcell_limiters.jl index 055e7ce24a4..554c4810178 100644 --- a/src/solvers/dgsem_tree/subcell_limiters.jl +++ b/src/solvers/dgsem_tree/subcell_limiters.jl @@ -16,18 +16,28 @@ end SubcellLimiterIDP(equations::AbstractEquations, basis; local_minmax_variables_cons = String[], positivity_variables_cons = String[], - positivity_correction_factor = 0.1) + positivity_variables_nonlinear = [], + positivity_correction_factor = 0.1, + max_iterations_newton = 10, + newton_tolerances = (1.0e-12, 1.0e-14), + gamma_constant_newton = 2 * ndims(equations)) Subcell invariant domain preserving (IDP) limiting used with [`VolumeIntegralSubcellLimiting`](@ref) including: - Local maximum/minimum Zalesak-type limiting for conservative variables (`local_minmax_variables_cons`) -- Positivity limiting for conservative variables (`positivity_variables_cons`) +- Positivity limiting for conservative variables (`positivity_variables_cons`) and non-linear variables +(`positivity_variables_nonlinear`) Conservative variables to be limited are passed as a vector of strings, e.g. `local_minmax_variables_cons = ["rho"]` -and `positivity_variables_cons = ["rho"]`. +and `positivity_variables_cons = ["rho"]`. For non-linear variables the specific functions are +passed in a vector, e.g. `positivity_variables_nonlinear = [pressure]`. The bounds are calculated using the low-order FV solution. The positivity limiter uses `positivity_correction_factor` such that `u^new >= positivity_correction_factor * u^FV`. +The limiting of non-linear variables uses a Newton-bisection method with a maximum of +`max_iterations_newton` iterations, relative and absolute tolerances of `newton_tolerances` +and a gamma constant of `gamma_constant_newton` (`gamma_constant_newton>=2*d`, +where `d = #dimensions`). !!! note This limiter and the correction callback [`SubcellLimiterIDPCorrection`](@ref) only work together. @@ -45,22 +55,32 @@ The bounds are calculated using the low-order FV solution. The positivity limite !!! warning "Experimental implementation" This is an experimental feature and may change in future releases. """ -struct SubcellLimiterIDP{RealT <: Real, Cache} <: AbstractSubcellLimiter +struct SubcellLimiterIDP{RealT <: Real, LimitingVariablesNonlinear, Cache} <: + AbstractSubcellLimiter local_minmax::Bool local_minmax_variables_cons::Vector{Int} # Local mininum/maximum principles for conservative variables positivity::Bool positivity_variables_cons::Vector{Int} # Positivity for conservative variables + positivity_variables_nonlinear::LimitingVariablesNonlinear # Positivity for nonlinear variables positivity_correction_factor::RealT cache::Cache + max_iterations_newton::Int + newton_tolerances::Tuple{RealT, RealT} # Relative and absolute tolerances for Newton's method + gamma_constant_newton::RealT # Constant for the subcell limiting of convex (nonlinear) constraints end # this method is used when the limiter is constructed as for shock-capturing volume integrals function SubcellLimiterIDP(equations::AbstractEquations, basis; local_minmax_variables_cons = String[], positivity_variables_cons = String[], - positivity_correction_factor = 0.1) + positivity_variables_nonlinear = [], + positivity_correction_factor = 0.1, + max_iterations_newton = 10, + newton_tolerances = (1.0e-12, 1.0e-14), + gamma_constant_newton = 2 * ndims(equations)) local_minmax = (length(local_minmax_variables_cons) > 0) - positivity = (length(positivity_variables_cons) > 0) + positivity = (length(positivity_variables_cons) + + length(positivity_variables_nonlinear) > 0) local_minmax_variables_cons_ = get_variable_index.(local_minmax_variables_cons, equations) @@ -80,13 +100,20 @@ function SubcellLimiterIDP(equations::AbstractEquations, basis; bound_keys = (bound_keys..., Symbol(string(v), "_min")) end end + for variable in positivity_variables_nonlinear + bound_keys = (bound_keys..., Symbol(string(variable), "_min")) + end cache = create_cache(SubcellLimiterIDP, equations, basis, bound_keys) SubcellLimiterIDP{typeof(positivity_correction_factor), + typeof(positivity_variables_nonlinear), typeof(cache)}(local_minmax, local_minmax_variables_cons_, positivity, positivity_variables_cons_, - positivity_correction_factor, cache) + positivity_variables_nonlinear, + positivity_correction_factor, cache, + max_iterations_newton, newton_tolerances, + gamma_constant_newton) end function Base.show(io::IO, limiter::SubcellLimiterIDP) @@ -124,7 +151,7 @@ function Base.show(io::IO, ::MIME"text/plain", limiter::SubcellLimiterIDP) ] end if positivity - string = "positivity for conservative variables $(limiter.positivity_variables_cons)" + string = "positivity for conservative variables $(limiter.positivity_variables_cons) and $(limiter.positivity_variables_nonlinear)" setup = [setup..., "" => string] setup = [ setup..., diff --git a/src/solvers/dgsem_tree/subcell_limiters_2d.jl b/src/solvers/dgsem_tree/subcell_limiters_2d.jl index bc69e55f264..2ff8c9e8d0c 100644 --- a/src/solvers/dgsem_tree/subcell_limiters_2d.jl +++ b/src/solvers/dgsem_tree/subcell_limiters_2d.jl @@ -229,6 +229,11 @@ end idp_positivity!(alpha, limiter, u, dt, semi, variable) end + # Nonlinear variables + for variable in limiter.positivity_variables_nonlinear + idp_positivity_nonlinear!(alpha, limiter, u, dt, semi, variable) + end + return nothing end @@ -292,4 +297,180 @@ end return nothing end + +@inline function idp_positivity_nonlinear!(alpha, limiter, u, dt, semi, variable) + mesh, equations, dg, cache = mesh_equations_solver_cache(semi) + (; positivity_correction_factor) = limiter + + (; variable_bounds) = limiter.cache.subcell_limiter_coefficients + var_min = variable_bounds[Symbol(string(variable), "_min")] + + @threaded for element in eachelement(dg, semi.cache) + for j in eachnode(dg), i in eachnode(dg) + # Compute bound + u_local = get_node_vars(u, equations, dg, i, j, element) + var = variable(u_local, equations) + if var < 0 + error("Safe $variable is not safe. element=$element, node: $i $j, value=$var") + end + var_min[i, j, element] = positivity_correction_factor * var + + # Perform Newton's bisection method to find new alpha + newton_loops_alpha!(alpha, var_min[i, j, element], u_local, i, j, element, + variable, initial_check_nonnegative, + final_check_nonnegative, + dt, mesh, equations, dg, cache, limiter) + end + end + + return nothing +end + +@inline function newton_loops_alpha!(alpha, bound, u, i, j, element, variable, + initial_check, final_check, dt, mesh, equations, + dg, cache, limiter) + (; inverse_weights) = dg.basis + (; antidiffusive_flux1_L, antidiffusive_flux2_L, antidiffusive_flux1_R, antidiffusive_flux2_R) = cache.antidiffusive_fluxes + inverse_jacobian = cache.elements.inverse_jacobian[element] + + (; gamma_constant_newton) = limiter + + # negative xi direction + antidiffusive_flux = gamma_constant_newton * inverse_jacobian * inverse_weights[i] * + get_node_vars(antidiffusive_flux1_R, equations, dg, i, j, + element) + newton_loop!(alpha, bound, u, i, j, element, variable, initial_check, final_check, + equations, dt, limiter, antidiffusive_flux) + + # positive xi direction + antidiffusive_flux = -gamma_constant_newton * inverse_jacobian * + inverse_weights[i] * + get_node_vars(antidiffusive_flux1_L, equations, dg, i + 1, j, + element) + newton_loop!(alpha, bound, u, i, j, element, variable, initial_check, final_check, + equations, dt, limiter, antidiffusive_flux) + + # negative eta direction + antidiffusive_flux = gamma_constant_newton * inverse_jacobian * inverse_weights[j] * + get_node_vars(antidiffusive_flux2_R, equations, dg, i, j, + element) + newton_loop!(alpha, bound, u, i, j, element, variable, initial_check, final_check, + equations, dt, limiter, antidiffusive_flux) + + # positive eta direction + antidiffusive_flux = -gamma_constant_newton * inverse_jacobian * + inverse_weights[j] * + get_node_vars(antidiffusive_flux2_L, equations, dg, i, j + 1, + element) + newton_loop!(alpha, bound, u, i, j, element, variable, initial_check, final_check, + equations, dt, limiter, antidiffusive_flux) + + return nothing +end + +@inline function newton_loop!(alpha, bound, u, i, j, element, variable, initial_check, + final_check, equations, dt, limiter, antidiffusive_flux) + newton_reltol, newton_abstol = limiter.newton_tolerances + + beta = 1 - alpha[i, j, element] + + beta_L = 0 # alpha = 1 + beta_R = beta # No higher beta (lower alpha) than the current one + + u_curr = u + beta * dt * antidiffusive_flux + + # If state is valid, perform initial check and return if correction is not needed + if is_valid_state(u_curr, equations) + as = goal_function(variable, bound, u_curr, equations) + + initial_check(bound, as, newton_abstol) && return nothing + end + + # Newton iterations + for iter in 1:(limiter.max_iterations_newton) + beta_old = beta + + # If the state is valid, evaluate d(goal)/d(beta) + if is_valid_state(u_curr, equations) + dSdbeta = dgoal_function(variable, u_curr, dt, antidiffusive_flux, + equations) + else # Otherwise, perform a bisection step + dSdbeta = 0 + end + + if dSdbeta != 0 + # Update beta with Newton's method + beta = beta - as / dSdbeta + end + + # Check bounds + if (beta < beta_L) || (beta > beta_R) || (dSdbeta == 0) || isnan(beta) + # Out of bounds, do a bisection step + beta = 0.5 * (beta_L + beta_R) + # Get new u + u_curr = u + beta * dt * antidiffusive_flux + + # If the state is invalid, finish bisection step without checking tolerance and iterate further + if !is_valid_state(u_curr, equations) + beta_R = beta + continue + end + + # Check new beta for condition and update bounds + as = goal_function(variable, bound, u_curr, equations) + if initial_check(bound, as, newton_abstol) + # New beta fulfills condition + beta_L = beta + else + # New beta does not fulfill condition + beta_R = beta + end + else + # Get new u + u_curr = u + beta * dt * antidiffusive_flux + + # If the state is invalid, redefine right bound without checking tolerance and iterate further + if !is_valid_state(u_curr, equations) + beta_R = beta + continue + end + + # Evaluate goal function + as = goal_function(variable, bound, u_curr, equations) + end + + # Check relative tolerance + if abs(beta_old - beta) <= newton_reltol + break + end + + # Check absolute tolerance + if final_check(bound, as, newton_abstol) + break + end + end + + new_alpha = 1 - beta + if alpha[i, j, element] > new_alpha + newton_abstol + error("Alpha is getting smaller. old: $(alpha[i, j, element]), new: $new_alpha") + else + alpha[i, j, element] = new_alpha + end + + return nothing +end + +# Initial checks +@inline initial_check_nonnegative(bound, goal, newton_abstol) = goal <= 0 + +# Goal and d(Goal)d(u) function +@inline goal_function(variable, bound, u, equations) = bound - variable(u, equations) +@inline function dgoal_function(variable, u, dt, antidiffusive_flux, equations) + -dot(variable(u, equations, True()), dt * antidiffusive_flux) +end + +# Final checks +@inline function final_check_nonnegative(bound, goal, newton_abstol) + (goal <= eps()) && (goal > -max(newton_abstol, abs(bound) * newton_abstol)) +end end # @muladd diff --git a/test/test_tree_2d_euler.jl b/test/test_tree_2d_euler.jl index 04a295537a3..dff78694a4a 100644 --- a/test/test_tree_2d_euler.jl +++ b/test/test_tree_2d_euler.jl @@ -581,6 +581,32 @@ end end end +@trixi_testset "elixir_euler_kelvin_helmholtz_instability_sc_subcell.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_euler_kelvin_helmholtz_instability_sc_subcell.jl"), + l2=[ + 0.05570313943221697, + 0.03298722119008495, + 0.05224470329798064, + 0.08011555326900793, + ], + linf=[ + 0.2409101307301882, + 0.16601896096921037, + 0.12356154796415487, + 0.2695166427148403, + ], + tspan=(0.0, 0.2)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 15000 + end +end + @trixi_testset "elixir_euler_colliding_flow.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_colliding_flow.jl"), l2=[ diff --git a/test/test_tree_2d_mhd.jl b/test/test_tree_2d_mhd.jl index 953c077c0a3..1f8458075aa 100644 --- a/test/test_tree_2d_mhd.jl +++ b/test/test_tree_2d_mhd.jl @@ -332,24 +332,28 @@ end @trixi_testset "elixir_mhd_shockcapturing_subcell.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_shockcapturing_subcell.jl"), - l2=[2.9974425783503109e-02, - 7.2849646345685956e-02, - 7.2488477174662239e-02, + l2=[ + 3.2064026219236076e-02, + 7.2461094392606618e-02, + 7.2380202888062711e-02, 0.0000000000000000e+00, - 1.2507971380965512e+00, - 1.8929505145499678e-02, - 1.2218606317164420e-02, + 8.6293936673145932e-01, + 8.4091669534557805e-03, + 5.2156364913231732e-03, 0.0000000000000000e+00, - 3.0154796910479838e-03], - linf=[3.2147382412340830e-01, - 1.3709471664007811e+00, - 1.3465154685288383e+00, + 2.0786952301129021e-04, + ], + linf=[ + 3.8778760255775635e-01, + 9.4666683953698927e-01, + 9.4618924645661928e-01, 0.0000000000000000e+00, - 1.6051257523415284e+01, - 3.0564266749926644e-01, - 2.3908016329805595e-01, + 1.0980297261521951e+01, + 1.0264404591009069e-01, + 1.0655686942176350e-01, 0.0000000000000000e+00, - 1.3711262178549158e-01], + 6.1013422157115546e-03, + ], tspan=(0.0, 0.003)) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) diff --git a/test/test_unit.jl b/test/test_unit.jl index 0c5f2bf0a4c..d0075e7ac93 100644 --- a/test/test_unit.jl +++ b/test/test_unit.jl @@ -414,7 +414,8 @@ end indicator_hg = IndicatorHennemannGassner(1.0, 0.0, true, "variable", "cache") @test_nowarn show(stdout, indicator_hg) - limiter_idp = SubcellLimiterIDP(true, [1], true, [1], 0.1, "cache") + limiter_idp = SubcellLimiterIDP(true, [1], true, [1], ["variable"], 0.1, "cache", 1, + (1.0, 1.0), 1.0) @test_nowarn show(stdout, limiter_idp) # TODO: TrixiShallowWater: move unit test From 849b09e6241dbc8732c5f8ecc49808490788401b Mon Sep 17 00:00:00 2001 From: bennibolm Date: Thu, 16 Nov 2023 12:17:44 +0100 Subject: [PATCH 02/25] Revise derivative function call; Add default derivative version --- src/equations/compressible_euler_2d.jl | 3 ++- src/equations/equations.jl | 6 ++++++ src/equations/ideal_glm_mhd_2d.jl | 3 ++- src/solvers/dgsem_tree/subcell_limiters_2d.jl | 2 +- 4 files changed, 11 insertions(+), 3 deletions(-) diff --git a/src/equations/compressible_euler_2d.jl b/src/equations/compressible_euler_2d.jl index f3d46da2451..ca10fbac665 100644 --- a/src/equations/compressible_euler_2d.jl +++ b/src/equations/compressible_euler_2d.jl @@ -1405,7 +1405,8 @@ end end # Transformation from conservative variables u to d(p)/d(u) -@inline function pressure(u, equations::CompressibleEulerEquations2D, derivative::True) +@inline function variable_derivative(::typeof(pressure), + u, equations::CompressibleEulerEquations2D) rho, rho_v1, rho_v2, rho_e = u v1 = rho_v1 / rho diff --git a/src/equations/equations.jl b/src/equations/equations.jl index ba2ad1cd1cd..fd8e858ae06 100644 --- a/src/equations/equations.jl +++ b/src/equations/equations.jl @@ -369,6 +369,12 @@ of the correct length `nvariables(equations)`. """ function energy_internal end +# Default implementation of derivation for `variable`. Used for subcell limiting. +# Implementing a derivative function for a specific function improves the performance. +@inline function variable_derivative(variable, u, equations) + return ForwardDiff.gradient(x -> variable(x, equations), u) +end + #################################################################################################### # Include files with actual implementations for different systems of equations. diff --git a/src/equations/ideal_glm_mhd_2d.jl b/src/equations/ideal_glm_mhd_2d.jl index dd80803e6a6..e486b5532b5 100644 --- a/src/equations/ideal_glm_mhd_2d.jl +++ b/src/equations/ideal_glm_mhd_2d.jl @@ -1063,7 +1063,8 @@ end end # Transformation from conservative variables u to d(p)/d(u) -@inline function pressure(u, equations::IdealGlmMhdEquations2D, derivative::True) +@inline function variable_derivative(::typeof(pressure), + u, equations::IdealGlmMhdEquations2D) rho, rho_v1, rho_v2, rho_v3, rho_e, B1, B2, B3, psi = u v1 = rho_v1 / rho diff --git a/src/solvers/dgsem_tree/subcell_limiters_2d.jl b/src/solvers/dgsem_tree/subcell_limiters_2d.jl index 2ff8c9e8d0c..e362d71cacc 100644 --- a/src/solvers/dgsem_tree/subcell_limiters_2d.jl +++ b/src/solvers/dgsem_tree/subcell_limiters_2d.jl @@ -466,7 +466,7 @@ end # Goal and d(Goal)d(u) function @inline goal_function(variable, bound, u, equations) = bound - variable(u, equations) @inline function dgoal_function(variable, u, dt, antidiffusive_flux, equations) - -dot(variable(u, equations, True()), dt * antidiffusive_flux) + -dot(variable_derivative(variable, u, equations), dt * antidiffusive_flux) end # Final checks From 2c75e383e768df60613a5b07605fcafebee24e82 Mon Sep 17 00:00:00 2001 From: bennibolm Date: Thu, 16 Nov 2023 16:32:59 +0100 Subject: [PATCH 03/25] Adapt test to actually test pos limiter for nonlinear variables --- test/test_tree_2d_euler.jl | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/test/test_tree_2d_euler.jl b/test/test_tree_2d_euler.jl index dff78694a4a..9315ae4d281 100644 --- a/test/test_tree_2d_euler.jl +++ b/test/test_tree_2d_euler.jl @@ -585,18 +585,18 @@ end @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_kelvin_helmholtz_instability_sc_subcell.jl"), l2=[ - 0.05570313943221697, - 0.03298722119008495, - 0.05224470329798064, - 0.08011555326900793, + 0.42185634563805724, + 0.1686471269704017, + 0.18240674916968103, + 0.17858250604280654, ], linf=[ - 0.2409101307301882, - 0.16601896096921037, - 0.12356154796415487, - 0.2695166427148403, + 1.7012978064377158, + 0.7149714986746726, + 0.5822547982757897, + 0.7300051017382696, ], - tspan=(0.0, 0.2)) + tspan=(0.0, 2.0)) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) let From bf69564e35851b037d86e921368bc1d0d6b18ecd Mon Sep 17 00:00:00 2001 From: bennibolm Date: Thu, 16 Nov 2023 17:07:51 +0100 Subject: [PATCH 04/25] Add unit test to test default implementation of variable_derivative --- test/test_unit.jl | 20 ++++++++++++++++++++ 1 file changed, 20 insertions(+) diff --git a/test/test_unit.jl b/test/test_unit.jl index d0075e7ac93..6888f62a615 100644 --- a/test/test_unit.jl +++ b/test/test_unit.jl @@ -1105,6 +1105,26 @@ end end end +@testset "Consistency check for variable_derivative routine" begin + # Set up conservative variables, equations + u = [ + 0.5011914484393387, + 0.8829127712445113, + 0.43024132987932817, + 0.7560616633050348, + ] + + equations = CompressibleEulerEquations2D(1.4) + + # Define wrapper function for pressure in order to call default implementation + function pressure_test(u, equations) + return pressure(u, equations) + end + + @test Trixi.variable_derivative(pressure_test, u, equations) ≈ + Trixi.variable_derivative(pressure, u, equations) +end + @testset "Equivalent Fluxes" begin # Set up equations and dummy conservative variables state # Burgers' Equation From 36b2481b2049313f24a10b8cd3231e3cd8ba73e3 Mon Sep 17 00:00:00 2001 From: bennibolm Date: Mon, 20 Nov 2023 13:14:42 +0100 Subject: [PATCH 05/25] Clean up comments and code --- src/solvers/dgsem_tree/subcell_limiters.jl | 18 ++++----- src/solvers/dgsem_tree/subcell_limiters_2d.jl | 38 +++++++++++++++---- 2 files changed, 39 insertions(+), 17 deletions(-) diff --git a/src/solvers/dgsem_tree/subcell_limiters.jl b/src/solvers/dgsem_tree/subcell_limiters.jl index 554c4810178..86d70e3853a 100644 --- a/src/solvers/dgsem_tree/subcell_limiters.jl +++ b/src/solvers/dgsem_tree/subcell_limiters.jl @@ -25,16 +25,16 @@ end Subcell invariant domain preserving (IDP) limiting used with [`VolumeIntegralSubcellLimiting`](@ref) including: - Local maximum/minimum Zalesak-type limiting for conservative variables (`local_minmax_variables_cons`) -- Positivity limiting for conservative variables (`positivity_variables_cons`) and non-linear variables +- Positivity limiting for conservative variables (`positivity_variables_cons`) and nonlinear variables (`positivity_variables_nonlinear`) Conservative variables to be limited are passed as a vector of strings, e.g. `local_minmax_variables_cons = ["rho"]` -and `positivity_variables_cons = ["rho"]`. For non-linear variables the specific functions are +and `positivity_variables_cons = ["rho"]`. For nonlinear variables the specific functions are passed in a vector, e.g. `positivity_variables_nonlinear = [pressure]`. The bounds are calculated using the low-order FV solution. The positivity limiter uses `positivity_correction_factor` such that `u^new >= positivity_correction_factor * u^FV`. -The limiting of non-linear variables uses a Newton-bisection method with a maximum of +The limiting of nonlinear variables uses a Newton-bisection method with a maximum of `max_iterations_newton` iterations, relative and absolute tolerances of `newton_tolerances` and a gamma constant of `gamma_constant_newton` (`gamma_constant_newton>=2*d`, where `d = #dimensions`). @@ -124,9 +124,9 @@ function Base.show(io::IO, limiter::SubcellLimiterIDP) if !(local_minmax || positivity) print(io, "No limiter selected => pure DG method") else - print(io, "limiter=(") - local_minmax && print(io, "min/max limiting, ") - positivity && print(io, "positivity") + print(io, "Limiter=(") + local_minmax && print(io, "Local min/max, ") + positivity && print(io, "Positivity, ") print(io, "), ") end print(io, "Local bounds with FV solution") @@ -147,15 +147,15 @@ function Base.show(io::IO, ::MIME"text/plain", limiter::SubcellLimiterIDP) if local_minmax setup = [ setup..., - "" => "local maximum/minimum bounds for conservative variables $(limiter.local_minmax_variables_cons)", + "" => "Local maximum/minimum limiting for conservative variables $(limiter.local_minmax_variables_cons)", ] end if positivity - string = "positivity for conservative variables $(limiter.positivity_variables_cons) and $(limiter.positivity_variables_nonlinear)" + string = "Positivity limiting for conservative variables $(limiter.positivity_variables_cons) and $(limiter.positivity_variables_nonlinear)" setup = [setup..., "" => string] setup = [ setup..., - "" => " positivity correction factor = $(limiter.positivity_correction_factor)", + "" => "- with positivity correction factor = $(limiter.positivity_correction_factor)", ] end setup = [ diff --git a/src/solvers/dgsem_tree/subcell_limiters_2d.jl b/src/solvers/dgsem_tree/subcell_limiters_2d.jl index e362d71cacc..1b6e8d318e7 100644 --- a/src/solvers/dgsem_tree/subcell_limiters_2d.jl +++ b/src/solvers/dgsem_tree/subcell_limiters_2d.jl @@ -5,6 +5,10 @@ @muladd begin #! format: noindent +############################################################################### +# IDP Limiting +############################################################################### + # this method is used when the limiter is constructed as for shock-capturing volume integrals function create_cache(limiter::Type{SubcellLimiterIDP}, equations::AbstractEquations{2}, basis::LobattoLegendreBasis, bound_keys) @@ -56,6 +60,9 @@ function (limiter::SubcellLimiterIDP)(u::AbstractArray{<:Any, 4}, semi, dg::DGSE return nothing end +############################################################################### +# Calculation of local bounds using low-order FV solution + @inline function calc_bounds_twosided!(var_min, var_max, variable, u, t, semi) mesh, equations, dg, cache = mesh_equations_solver_cache(semi) # Calc bounds inside elements @@ -154,6 +161,9 @@ end return nothing end +############################################################################### +# Local minimum/maximum limiting + @inline function idp_local_minmax!(alpha, limiter, u, t, dt, semi) for variable in limiter.local_minmax_variables_cons idp_local_minmax!(alpha, limiter, u, t, dt, semi, variable) @@ -223,20 +233,32 @@ end return nothing end +############################################################################### +# Global positivity limiting + @inline function idp_positivity!(alpha, limiter, u, dt, semi) # Conservative variables for variable in limiter.positivity_variables_cons - idp_positivity!(alpha, limiter, u, dt, semi, variable) + @trixi_timeit timer() "conservative variables" idp_positivity!(alpha, limiter, + u, dt, semi, + variable) end # Nonlinear variables for variable in limiter.positivity_variables_nonlinear - idp_positivity_nonlinear!(alpha, limiter, u, dt, semi, variable) + @trixi_timeit timer() "nonlinear variables" idp_positivity_nonlinear!(alpha, + limiter, + u, dt, + semi, + variable) end return nothing end +############################################################################### +# Global positivity limiting of conservative variables + @inline function idp_positivity!(alpha, limiter, u, dt, semi, variable) mesh, equations, dg, cache = mesh_equations_solver_cache(semi) (; antidiffusive_flux1_L, antidiffusive_flux2_L, antidiffusive_flux1_R, antidiffusive_flux2_R) = cache.antidiffusive_fluxes @@ -299,13 +321,14 @@ end end @inline function idp_positivity_nonlinear!(alpha, limiter, u, dt, semi, variable) - mesh, equations, dg, cache = mesh_equations_solver_cache(semi) + _, equations, dg, cache = mesh_equations_solver_cache(semi) (; positivity_correction_factor) = limiter (; variable_bounds) = limiter.cache.subcell_limiter_coefficients var_min = variable_bounds[Symbol(string(variable), "_min")] @threaded for element in eachelement(dg, semi.cache) + inverse_jacobian = cache.elements.inverse_jacobian[element] for j in eachnode(dg), i in eachnode(dg) # Compute bound u_local = get_node_vars(u, equations, dg, i, j, element) @@ -318,8 +341,8 @@ end # Perform Newton's bisection method to find new alpha newton_loops_alpha!(alpha, var_min[i, j, element], u_local, i, j, element, variable, initial_check_nonnegative, - final_check_nonnegative, - dt, mesh, equations, dg, cache, limiter) + final_check_nonnegative, inverse_jacobian, + dt, equations, dg, cache, limiter) end end @@ -327,11 +350,10 @@ end end @inline function newton_loops_alpha!(alpha, bound, u, i, j, element, variable, - initial_check, final_check, dt, mesh, equations, - dg, cache, limiter) + initial_check, final_check, inverse_jacobian, dt, + equations, dg, cache, limiter) (; inverse_weights) = dg.basis (; antidiffusive_flux1_L, antidiffusive_flux2_L, antidiffusive_flux1_R, antidiffusive_flux2_R) = cache.antidiffusive_fluxes - inverse_jacobian = cache.elements.inverse_jacobian[element] (; gamma_constant_newton) = limiter From 2702bb8ae9422701bf7279d46f25fb1a9bdeb0fa Mon Sep 17 00:00:00 2001 From: bennibolm Date: Thu, 7 Dec 2023 16:48:16 +0100 Subject: [PATCH 06/25] Rename Newton-bisection variables --- src/solvers/dgsem_tree/subcell_limiters_2d.jl | 24 +++++++++---------- 1 file changed, 12 insertions(+), 12 deletions(-) diff --git a/src/solvers/dgsem_tree/subcell_limiters_2d.jl b/src/solvers/dgsem_tree/subcell_limiters_2d.jl index fee22e188ad..788d32344a7 100644 --- a/src/solvers/dgsem_tree/subcell_limiters_2d.jl +++ b/src/solvers/dgsem_tree/subcell_limiters_2d.jl @@ -402,9 +402,9 @@ end # If state is valid, perform initial check and return if correction is not needed if is_valid_state(u_curr, equations) - as = goal_function(variable, bound, u_curr, equations) + goal = goal_function(variable, bound, u_curr, equations) - initial_check(bound, as, newton_abstol) && return nothing + initial_check(bound, goal, newton_abstol) && return nothing end # Newton iterations @@ -413,19 +413,19 @@ end # If the state is valid, evaluate d(goal)/d(beta) if is_valid_state(u_curr, equations) - dSdbeta = dgoal_function(variable, u_curr, dt, antidiffusive_flux, - equations) + dgoal_dbeta = dgoal_function(variable, u_curr, dt, antidiffusive_flux, + equations) else # Otherwise, perform a bisection step - dSdbeta = 0 + dgoal_dbeta = 0 end - if dSdbeta != 0 + if dgoal_dbeta != 0 # Update beta with Newton's method - beta = beta - as / dSdbeta + beta = beta - goal / dgoal_dbeta end # Check bounds - if (beta < beta_L) || (beta > beta_R) || (dSdbeta == 0) || isnan(beta) + if (beta < beta_L) || (beta > beta_R) || (dgoal_dbeta == 0) || isnan(beta) # Out of bounds, do a bisection step beta = 0.5 * (beta_L + beta_R) # Get new u @@ -438,8 +438,8 @@ end end # Check new beta for condition and update bounds - as = goal_function(variable, bound, u_curr, equations) - if initial_check(bound, as, newton_abstol) + goal = goal_function(variable, bound, u_curr, equations) + if initial_check(bound, goal, newton_abstol) # New beta fulfills condition beta_L = beta else @@ -457,7 +457,7 @@ end end # Evaluate goal function - as = goal_function(variable, bound, u_curr, equations) + goal = goal_function(variable, bound, u_curr, equations) end # Check relative tolerance @@ -466,7 +466,7 @@ end end # Check absolute tolerance - if final_check(bound, as, newton_abstol) + if final_check(bound, goal, newton_abstol) break end end From 5ec5adee3d248bb255868a32b9790ea717a42c2f Mon Sep 17 00:00:00 2001 From: bennibolm Date: Wed, 13 Dec 2023 12:59:26 +0100 Subject: [PATCH 07/25] Implement suggestions --- src/solvers/dgsem_tree/subcell_limiters.jl | 4 ++-- src/solvers/dgsem_tree/subcell_limiters_2d.jl | 11 +++++++---- 2 files changed, 9 insertions(+), 6 deletions(-) diff --git a/src/solvers/dgsem_tree/subcell_limiters.jl b/src/solvers/dgsem_tree/subcell_limiters.jl index 86d70e3853a..fb52272e2e8 100644 --- a/src/solvers/dgsem_tree/subcell_limiters.jl +++ b/src/solvers/dgsem_tree/subcell_limiters.jl @@ -36,8 +36,8 @@ The bounds are calculated using the low-order FV solution. The positivity limite `positivity_correction_factor` such that `u^new >= positivity_correction_factor * u^FV`. The limiting of nonlinear variables uses a Newton-bisection method with a maximum of `max_iterations_newton` iterations, relative and absolute tolerances of `newton_tolerances` -and a gamma constant of `gamma_constant_newton` (`gamma_constant_newton>=2*d`, -where `d = #dimensions`). +and a provisional update constant `gamma_constant_newton` (`gamma_constant_newton>=2*d`, +where `d = #dimensions`). See equation (20) of Pazner (2020) and equation (30) of Rueda-Ramírez et al. (2022). !!! note This limiter and the correction callback [`SubcellLimiterIDPCorrection`](@ref) only work together. diff --git a/src/solvers/dgsem_tree/subcell_limiters_2d.jl b/src/solvers/dgsem_tree/subcell_limiters_2d.jl index 788d32344a7..98d2a3296cd 100644 --- a/src/solvers/dgsem_tree/subcell_limiters_2d.jl +++ b/src/solvers/dgsem_tree/subcell_limiters_2d.jl @@ -238,9 +238,12 @@ end @inline function idp_positivity!(alpha, limiter, u, dt, semi) # Conservative variables for variable in limiter.positivity_variables_cons - @trixi_timeit timer() "conservative variables" idp_positivity!(alpha, limiter, - u, dt, semi, - variable) + @trixi_timeit timer() "conservative variables" idp_positivity_conservative!(alpha, + limiter, + u, + dt, + semi, + variable) end # Nonlinear variables @@ -258,7 +261,7 @@ end ############################################################################### # Global positivity limiting of conservative variables -@inline function idp_positivity!(alpha, limiter, u, dt, semi, variable) +@inline function idp_positivity_conservative!(alpha, limiter, u, dt, semi, variable) mesh, equations, dg, cache = mesh_equations_solver_cache(semi) (; antidiffusive_flux1_L, antidiffusive_flux2_L, antidiffusive_flux1_R, antidiffusive_flux2_R) = cache.antidiffusive_fluxes (; inverse_weights) = dg.basis From ccdc34edf591ac5bab10fb4b2011c86da6b3b6ac Mon Sep 17 00:00:00 2001 From: bennibolm Date: Fri, 22 Dec 2023 12:27:40 +0100 Subject: [PATCH 08/25] Relocate functions --- src/equations/compressible_euler_2d.jl | 41 +++++++++++------------ src/equations/ideal_glm_mhd_2d.jl | 45 +++++++++++++------------- 2 files changed, 44 insertions(+), 42 deletions(-) diff --git a/src/equations/compressible_euler_2d.jl b/src/equations/compressible_euler_2d.jl index ca10fbac665..47916283a17 100644 --- a/src/equations/compressible_euler_2d.jl +++ b/src/equations/compressible_euler_2d.jl @@ -1404,18 +1404,6 @@ end return SVector(w1, w2, w3, w4) end -# Transformation from conservative variables u to d(p)/d(u) -@inline function variable_derivative(::typeof(pressure), - u, equations::CompressibleEulerEquations2D) - rho, rho_v1, rho_v2, rho_e = u - - v1 = rho_v1 / rho - v2 = rho_v2 / rho - v_square = v1^2 + v2^2 - - return (equations.gamma - 1.0) * SVector(0.5 * v_square, -v1, -v2, 1.0) -end - @inline function entropy2cons(w, equations::CompressibleEulerEquations2D) # See Hughes, Franca, Mallet (1986) A new finite element formulation for CFD # [DOI: 10.1016/0045-7825(86)90127-1](https://doi.org/10.1016/0045-7825(86)90127-1) @@ -1440,14 +1428,6 @@ end return SVector(rho, rho_v1, rho_v2, rho_e) end -@inline function is_valid_state(cons, equations::CompressibleEulerEquations2D) - p = pressure(cons, equations) - if cons[1] <= 0.0 || p <= 0.0 - return false - end - return true -end - # Convert primitive to conservative variables @inline function prim2cons(prim, equations::CompressibleEulerEquations2D) rho, v1, v2, p = prim @@ -1468,6 +1448,18 @@ end return p end +# Transformation from conservative variables u to d(p)/d(u) +@inline function variable_derivative(::typeof(pressure), + u, equations::CompressibleEulerEquations2D) + rho, rho_v1, rho_v2, rho_e = u + + v1 = rho_v1 / rho + v2 = rho_v2 / rho + v_square = v1^2 + v2^2 + + return (equations.gamma - 1.0) * SVector(0.5 * v_square, -v1, -v2, 1.0) +end + @inline function density_pressure(u, equations::CompressibleEulerEquations2D) rho, rho_v1, rho_v2, rho_e = u rho_times_p = (equations.gamma - 1) * (rho * rho_e - 0.5 * (rho_v1^2 + rho_v2^2)) @@ -1535,4 +1527,13 @@ end @inline function energy_internal(cons, equations::CompressibleEulerEquations2D) return energy_total(cons, equations) - energy_kinetic(cons, equations) end + +# State validation for subcell limiting using Newton-bisection method +@inline function is_valid_state(cons, equations::CompressibleEulerEquations2D) + p = pressure(cons, equations) + if cons[1] <= 0.0 || p <= 0.0 + return false + end + return true +end end # @muladd diff --git a/src/equations/ideal_glm_mhd_2d.jl b/src/equations/ideal_glm_mhd_2d.jl index e486b5532b5..b1371ba38b3 100644 --- a/src/equations/ideal_glm_mhd_2d.jl +++ b/src/equations/ideal_glm_mhd_2d.jl @@ -1062,20 +1062,6 @@ end return SVector(w1, w2, w3, w4, w5, w6, w7, w8, w9) end -# Transformation from conservative variables u to d(p)/d(u) -@inline function variable_derivative(::typeof(pressure), - u, equations::IdealGlmMhdEquations2D) - rho, rho_v1, rho_v2, rho_v3, rho_e, B1, B2, B3, psi = u - - v1 = rho_v1 / rho - v2 = rho_v2 / rho - v3 = rho_v3 / rho - v_square = v1^2 + v2^2 + v3^2 - - return (equations.gamma - 1.0) * - SVector(0.5 * v_square, -v1, -v2, -v3, 1.0, -B1, -B2, -B3, -psi) -end - # Convert entropy variables to conservative variables @inline function entropy2cons(w, equations::IdealGlmMhdEquations2D) w1, w2, w3, w4, w5, w6, w7, w8, w9 = w @@ -1103,14 +1089,6 @@ end return prim2cons(SVector(rho, v1, v2, v3, p, B1, B2, B3, psi), equations) end -@inline function is_valid_state(cons, equations::IdealGlmMhdEquations2D) - p = pressure(cons, equations) - if cons[1] <= 0.0 || p <= 0.0 - return false - end - return true -end - # Convert primitive to conservative variables @inline function prim2cons(prim, equations::IdealGlmMhdEquations2D) rho, v1, v2, v3, p, B1, B2, B3, psi = prim @@ -1140,6 +1118,20 @@ end return p end +# Transformation from conservative variables u to d(p)/d(u) +@inline function variable_derivative(::typeof(pressure), + u, equations::IdealGlmMhdEquations2D) + rho, rho_v1, rho_v2, rho_v3, rho_e, B1, B2, B3, psi = u + + v1 = rho_v1 / rho + v2 = rho_v2 / rho + v3 = rho_v3 / rho + v_square = v1^2 + v2^2 + v3^2 + + return (equations.gamma - 1.0) * + SVector(0.5 * v_square, -v1, -v2, -v3, 1.0, -B1, -B2, -B3, -psi) +end + @inline function density_pressure(u, equations::IdealGlmMhdEquations2D) rho, rho_v1, rho_v2, rho_v3, rho_e, B1, B2, B3, psi = u p = (equations.gamma - 1) * (rho_e - 0.5 * (rho_v1^2 + rho_v2^2 + rho_v3^2) / rho @@ -1406,6 +1398,15 @@ end cons[9]^2 / 2) end +# State validation for subcell limiting using Newton-bisection method +@inline function is_valid_state(cons, equations::IdealGlmMhdEquations2D) + p = pressure(cons, equations) + if cons[1] <= 0.0 || p <= 0.0 + return false + end + return true +end + # Calculate the cross helicity (\vec{v}⋅\vec{B}) for a conservative state `cons' @inline function cross_helicity(cons, ::IdealGlmMhdEquations2D) return (cons[2] * cons[6] + cons[3] * cons[7] + cons[4] * cons[8]) / cons[1] From 91e5239375b164cf61e1a83c3ee192328e91e8fa Mon Sep 17 00:00:00 2001 From: bennibolm Date: Fri, 22 Dec 2023 13:51:23 +0100 Subject: [PATCH 09/25] Add entropy limiters --- ...euler_blast_wave_sc_subcell_nonperiodic.jl | 3 +- ...lixir_euler_sedov_blast_wave_sc_subcell.jl | 3 +- src/callbacks_stage/subcell_bounds_check.jl | 18 +- .../subcell_bounds_check_2d.jl | 32 +++- src/equations/compressible_euler_2d.jl | 52 ++++++ src/solvers/dgsem_tree/subcell_limiters.jl | 43 ++++- src/solvers/dgsem_tree/subcell_limiters_2d.jl | 169 +++++++++++++++++- test/test_tree_2d_euler.jl | 32 ++-- test/test_unit.jl | 4 +- 9 files changed, 325 insertions(+), 31 deletions(-) diff --git a/examples/tree_2d_dgsem/elixir_euler_blast_wave_sc_subcell_nonperiodic.jl b/examples/tree_2d_dgsem/elixir_euler_blast_wave_sc_subcell_nonperiodic.jl index 209aa2ae352..6ff381b896c 100644 --- a/examples/tree_2d_dgsem/elixir_euler_blast_wave_sc_subcell_nonperiodic.jl +++ b/examples/tree_2d_dgsem/elixir_euler_blast_wave_sc_subcell_nonperiodic.jl @@ -41,7 +41,8 @@ surface_flux = flux_lax_friedrichs volume_flux = flux_ranocha basis = LobattoLegendreBasis(3) limiter_idp = SubcellLimiterIDP(equations, basis; - local_minmax_variables_cons = ["rho"]) + local_minmax_variables_cons = ["rho"], + math_entropy = true) volume_integral = VolumeIntegralSubcellLimiting(limiter_idp; volume_flux_dg = volume_flux, volume_flux_fv = surface_flux) diff --git a/examples/tree_2d_dgsem/elixir_euler_sedov_blast_wave_sc_subcell.jl b/examples/tree_2d_dgsem/elixir_euler_sedov_blast_wave_sc_subcell.jl index c1ba3d96962..64cae963158 100644 --- a/examples/tree_2d_dgsem/elixir_euler_sedov_blast_wave_sc_subcell.jl +++ b/examples/tree_2d_dgsem/elixir_euler_sedov_blast_wave_sc_subcell.jl @@ -42,7 +42,8 @@ surface_flux = flux_lax_friedrichs volume_flux = flux_chandrashekar basis = LobattoLegendreBasis(3) limiter_idp = SubcellLimiterIDP(equations, basis; - local_minmax_variables_cons = ["rho"]) + local_minmax_variables_cons = ["rho"], + spec_entropy = true) volume_integral = VolumeIntegralSubcellLimiting(limiter_idp; volume_flux_dg = volume_flux, volume_flux_fv = surface_flux) diff --git a/src/callbacks_stage/subcell_bounds_check.jl b/src/callbacks_stage/subcell_bounds_check.jl index 25662f85802..04fe084effd 100644 --- a/src/callbacks_stage/subcell_bounds_check.jl +++ b/src/callbacks_stage/subcell_bounds_check.jl @@ -77,7 +77,7 @@ function init_callback(callback::BoundsCheckCallback, semi, limiter::SubcellLimi return nothing end - (; local_minmax, positivity) = limiter + (; local_minmax, positivity, spec_entropy, math_entropy) = limiter (; output_directory) = callback variables = varnames(cons2cons, semi.equations) @@ -90,6 +90,12 @@ function init_callback(callback::BoundsCheckCallback, semi, limiter::SubcellLimi print(f, ", " * variable_string * "_min, " * variable_string * "_max") end end + if spec_entropy + print(f, ", specEntr_min") + end + if math_entropy + print(f, ", mathEntr_max") + end if positivity for v in limiter.positivity_variables_cons if v in limiter.local_minmax_variables_cons @@ -120,7 +126,7 @@ end @inline function finalize_callback(callback::BoundsCheckCallback, semi, limiter::SubcellLimiterIDP) - (; local_minmax, positivity) = limiter + (; local_minmax, positivity, spec_entropy, math_entropy) = limiter (; idp_bounds_delta) = limiter.cache variables = varnames(cons2cons, semi.equations) @@ -135,6 +141,14 @@ end println("-upper bound: ", idp_bounds_delta[Symbol(v_string, "_max")][2]) end end + if spec_entropy + println("spec. entropy:\n- lower bound: ", + idp_bounds_delta[:spec_entropy_min][2]) + end + if math_entropy + println("math. entropy:\n- upper bound: ", + idp_bounds_delta[:math_entropy_max][2]) + end if positivity for v in limiter.positivity_variables_cons if v in limiter.local_minmax_variables_cons diff --git a/src/callbacks_stage/subcell_bounds_check_2d.jl b/src/callbacks_stage/subcell_bounds_check_2d.jl index 688d3281b9e..06e5663e120 100644 --- a/src/callbacks_stage/subcell_bounds_check_2d.jl +++ b/src/callbacks_stage/subcell_bounds_check_2d.jl @@ -8,7 +8,7 @@ @inline function check_bounds(u, mesh::AbstractMesh{2}, equations, solver, cache, limiter::SubcellLimiterIDP, time, iter, output_directory, save_errors) - (; local_minmax, positivity) = solver.volume_integral.limiter + (; local_minmax, positivity, spec_entropy, math_entropy) = solver.volume_integral.limiter (; variable_bounds) = limiter.cache.subcell_limiter_coefficients (; idp_bounds_delta) = limiter.cache @@ -32,6 +32,30 @@ deviation_max[2] = max(deviation_max[2], deviation_max[1]) end end + if spec_entropy + key = :spec_entropy_min + deviation = idp_bounds_delta[key] + for element in eachelement(solver, cache), j in eachnode(solver), + i in eachnode(solver) + + s = entropy_spec(get_node_vars(u, equations, solver, i, j, element), + equations) + deviation[1] = max(deviation[1], variable_bounds[key][i, j, element] - s) + end + deviation[2] = max(deviation[2], deviation[1]) + end + if math_entropy + key = :math_entropy_max + deviation = idp_bounds_delta[key] + for element in eachelement(solver, cache), j in eachnode(solver), + i in eachnode(solver) + + s = entropy_math(get_node_vars(u, equations, solver, i, j, element), + equations) + deviation[1] = max(deviation[1], s - variable_bounds[key][i, j, element]) + end + deviation[2] = max(deviation[2], deviation[1]) + end if positivity for v in limiter.positivity_variables_cons if v in limiter.local_minmax_variables_cons @@ -73,6 +97,12 @@ idp_bounds_delta[Symbol(v_string, "_max")][1]) end end + if spec_entropy + print(f, ", ", idp_bounds_delta[:spec_entropy_min][1]) + end + if math_entropy + print(f, ", ", idp_bounds_delta[:math_entropy_max][1]) + end if positivity for v in limiter.positivity_variables_cons if v in limiter.local_minmax_variables_cons diff --git a/src/equations/compressible_euler_2d.jl b/src/equations/compressible_euler_2d.jl index 47916283a17..d660f6f0db8 100644 --- a/src/equations/compressible_euler_2d.jl +++ b/src/equations/compressible_euler_2d.jl @@ -1404,6 +1404,33 @@ end return SVector(w1, w2, w3, w4) end +# Transformation from conservative variables u to entropy vector dSdu, +# S = -rho*s/(gamma-1), s=ln(p)-gamma*ln(rho) +@inline function cons2entropy_spec(u, equations::CompressibleEulerEquations2D) + rho, rho_v1, rho_v2, rho_e = u + + v1 = rho_v1 / rho + v2 = rho_v2 / rho + v_square = v1^2 + v2^2 + inv_rho_gammap1 = (1 / rho)^(equations.gamma + 1.0) + + # The derivative vector for the modified specific entropy of Guermond et al. + w1 = inv_rho_gammap1 * + (0.5 * rho * (equations.gamma + 1.0) * v_square - equations.gamma * rho_e) + w2 = -rho_v1 * inv_rho_gammap1 + w3 = -rho_v2 * inv_rho_gammap1 + w4 = (1 / rho)^equations.gamma + + # The derivative vector for other specific entropy + # sp = 1.0/(gammam1 * (rho_e - 0.5 * rho * v_square) + # w1 = gammam1 * 0.5 * v_square * sp - gamma / rho + # w2 = -gammam1 * v1 * sp + # w3 = -gammam1 * v2 * sp + # w4 = gammam1 * sp + + return SVector(w1, w2, w3, w4) +end + @inline function entropy2cons(w, equations::CompressibleEulerEquations2D) # See Hughes, Franca, Mallet (1986) A new finite element formulation for CFD # [DOI: 10.1016/0045-7825(86)90127-1](https://doi.org/10.1016/0045-7825(86)90127-1) @@ -1509,6 +1536,31 @@ end return S end +# Transformation from conservative variables u to d(s)/d(u) +@inline function variable_derivative(::typeof(entropy_math), + u, equations::CompressibleEulerEquations2D) + return cons2entropy(u, equations) +end + +# Calculate specific entropy for conservative variable u +@inline function entropy_spec(u, equations::CompressibleEulerEquations2D) + rho, rho_v1, rho_v2, rho_e = u + + # Modified specific entropy from Guermond et al. (2019) + s = (rho_e - 0.5 * (rho_v1^2 + rho_v2^2) / rho) * (1 / rho)^equations.gamma + + # Other specific entropy + # rho_sp = rho/((equations.gamma - 1.0) * (rho_e - 0.5 * rho * v_square)) + # s = log(p) - (equaions.gamma + 1) * log(rho) + return s +end + +# Transformation from conservative variables u to d(s)/d(u) +@inline function variable_derivative(::typeof(entropy_spec), + u, equations::CompressibleEulerEquations2D) + return cons2entropy_spec(u, equations) +end + # Default entropy is the mathematical entropy @inline function entropy(cons, equations::CompressibleEulerEquations2D) entropy_math(cons, equations) diff --git a/src/solvers/dgsem_tree/subcell_limiters.jl b/src/solvers/dgsem_tree/subcell_limiters.jl index fb52272e2e8..db443853845 100644 --- a/src/solvers/dgsem_tree/subcell_limiters.jl +++ b/src/solvers/dgsem_tree/subcell_limiters.jl @@ -18,6 +18,8 @@ end positivity_variables_cons = String[], positivity_variables_nonlinear = [], positivity_correction_factor = 0.1, + spec_entropy = false, + math_entropy = false, max_iterations_newton = 10, newton_tolerances = (1.0e-12, 1.0e-14), gamma_constant_newton = 2 * ndims(equations)) @@ -27,6 +29,7 @@ including: - Local maximum/minimum Zalesak-type limiting for conservative variables (`local_minmax_variables_cons`) - Positivity limiting for conservative variables (`positivity_variables_cons`) and nonlinear variables (`positivity_variables_nonlinear`) +- One-sided limiting for specific and mathematical entropy (`spec_entropy`, `math_entropy`) Conservative variables to be limited are passed as a vector of strings, e.g. `local_minmax_variables_cons = ["rho"]` and `positivity_variables_cons = ["rho"]`. For nonlinear variables the specific functions are @@ -63,6 +66,8 @@ struct SubcellLimiterIDP{RealT <: Real, LimitingVariablesNonlinear, Cache} <: positivity_variables_cons::Vector{Int} # Positivity for conservative variables positivity_variables_nonlinear::LimitingVariablesNonlinear # Positivity for nonlinear variables positivity_correction_factor::RealT + spec_entropy::Bool + math_entropy::Bool cache::Cache max_iterations_newton::Int newton_tolerances::Tuple{RealT, RealT} # Relative and absolute tolerances for Newton's method @@ -75,12 +80,17 @@ function SubcellLimiterIDP(equations::AbstractEquations, basis; positivity_variables_cons = String[], positivity_variables_nonlinear = [], positivity_correction_factor = 0.1, + spec_entropy = false, + math_entropy = false, max_iterations_newton = 10, newton_tolerances = (1.0e-12, 1.0e-14), gamma_constant_newton = 2 * ndims(equations)) local_minmax = (length(local_minmax_variables_cons) > 0) positivity = (length(positivity_variables_cons) + length(positivity_variables_nonlinear) > 0) + if math_entropy && spec_entropy + error("Only one of the two can be selected: math_entropy/spec_entropy") + end local_minmax_variables_cons_ = get_variable_index.(local_minmax_variables_cons, equations) @@ -95,6 +105,12 @@ function SubcellLimiterIDP(equations::AbstractEquations, basis; Symbol(v_string, "_max")) end end + if spec_entropy + bound_keys = (bound_keys..., :spec_entropy_min) + end + if math_entropy + bound_keys = (bound_keys..., :math_entropy_max) + end for v in positivity_variables_cons_ if !(v in local_minmax_variables_cons_) bound_keys = (bound_keys..., Symbol(string(v), "_min")) @@ -111,22 +127,26 @@ function SubcellLimiterIDP(equations::AbstractEquations, basis; typeof(cache)}(local_minmax, local_minmax_variables_cons_, positivity, positivity_variables_cons_, positivity_variables_nonlinear, - positivity_correction_factor, cache, + positivity_correction_factor, + spec_entropy, math_entropy, + cache, max_iterations_newton, newton_tolerances, gamma_constant_newton) end function Base.show(io::IO, limiter::SubcellLimiterIDP) @nospecialize limiter # reduce precompilation time - (; local_minmax, positivity) = limiter + (; local_minmax, positivity, spec_entropy, math_entropy) = limiter print(io, "SubcellLimiterIDP(") - if !(local_minmax || positivity) + if !(local_minmax || positivity || spec_entropy || math_entropy) print(io, "No limiter selected => pure DG method") else print(io, "Limiter=(") local_minmax && print(io, "Local min/max, ") positivity && print(io, "Positivity, ") + spec_entropy && print(io, "Specific entropy, ") + math_entropy && print(io, "Mathematical entropy, ") print(io, "), ") end print(io, "Local bounds with FV solution") @@ -135,15 +155,15 @@ end function Base.show(io::IO, ::MIME"text/plain", limiter::SubcellLimiterIDP) @nospecialize limiter # reduce precompilation time - (; local_minmax, positivity) = limiter + (; local_minmax, positivity, spec_entropy, math_entropy) = limiter if get(io, :compact, false) show(io, limiter) else - if !(local_minmax || positivity) - setup = ["limiter" => "No limiter selected => pure DG method"] + if !(local_minmax || positivity || spec_entropy || math_entropy) + setup = ["Limiter" => "No limiter selected => pure DG method"] else - setup = ["limiter" => ""] + setup = ["Limiter" => ""] if local_minmax setup = [ setup..., @@ -158,6 +178,15 @@ function Base.show(io::IO, ::MIME"text/plain", limiter::SubcellLimiterIDP) "" => "- with positivity correction factor = $(limiter.positivity_correction_factor)", ] end + if spec_entropy + setup = [setup..., "" => "Local minimum limiting for specific entropy"] + end + if math_entropy + setup = [ + setup..., + "" => "Local maximum limiting for mathematical entropy", + ] + end setup = [ setup..., "Local bounds" => "FV solution", diff --git a/src/solvers/dgsem_tree/subcell_limiters_2d.jl b/src/solvers/dgsem_tree/subcell_limiters_2d.jl index 98d2a3296cd..edc2a2182ad 100644 --- a/src/solvers/dgsem_tree/subcell_limiters_2d.jl +++ b/src/solvers/dgsem_tree/subcell_limiters_2d.jl @@ -40,6 +40,14 @@ function (limiter::SubcellLimiterIDP)(u::AbstractArray{<:Any, 4}, semi, dg::DGSE if limiter.positivity @trixi_timeit timer() "positivity" idp_positivity!(alpha, limiter, u, dt, semi) end + if limiter.spec_entropy + @trixi_timeit timer() "spec_entropy" idp_spec_entropy!(alpha, limiter, u, t, dt, + semi) + end + if limiter.math_entropy + @trixi_timeit timer() "math_entropy" idp_math_entropy!(alpha, limiter, u, t, dt, + semi) + end # Calculate alpha1 and alpha2 @unpack alpha1, alpha2 = limiter.cache.subcell_limiter_coefficients @@ -160,6 +168,105 @@ end return nothing end +@inline function calc_bounds_onesided!(var_minmax, minmax, typeminmax, variable, u, t, + semi) + mesh, equations, dg, cache = mesh_equations_solver_cache(semi) + # Calc bounds inside elements + @threaded for element in eachelement(dg, cache) + for j in eachnode(dg), i in eachnode(dg) + var_minmax[i, j, element] = typeminmax(eltype(var_minmax)) + end + + # Calculate bounds at Gauss-Lobatto nodes using u + for j in eachnode(dg), i in eachnode(dg) + var = variable(get_node_vars(u, equations, dg, i, j, element), equations) + var_minmax[i, j, element] = minmax(var_minmax[i, j, element], var) + + if i > 1 + var_minmax[i - 1, j, element] = minmax(var_minmax[i - 1, j, element], + var) + end + if i < nnodes(dg) + var_minmax[i + 1, j, element] = minmax(var_minmax[i + 1, j, element], + var) + end + if j > 1 + var_minmax[i, j - 1, element] = minmax(var_minmax[i, j - 1, element], + var) + end + if j < nnodes(dg) + var_minmax[i, j + 1, element] = minmax(var_minmax[i, j + 1, element], + var) + end + end + end + + # Values at element boundary + calc_bounds_onesided_interface!(var_minmax, minmax, variable, u, t, semi, mesh) +end + +@inline function calc_bounds_onesided_interface!(var_minmax, minmax, variable, u, t, + semi, mesh::TreeMesh2D) + _, equations, dg, cache = mesh_equations_solver_cache(semi) + (; boundary_conditions) = semi + # Calc bounds at interfaces and periodic boundaries + for interface in eachinterface(dg, cache) + # Get neighboring element ids + left = cache.interfaces.neighbor_ids[1, interface] + right = cache.interfaces.neighbor_ids[2, interface] + + orientation = cache.interfaces.orientations[interface] + + for i in eachnode(dg) + index_left = (nnodes(dg), i) + index_right = (1, i) + if orientation == 2 + index_left = reverse(index_left) + index_right = reverse(index_right) + end + var_left = variable(get_node_vars(u, equations, dg, index_left..., left), + equations) + var_right = variable(get_node_vars(u, equations, dg, index_right..., right), + equations) + + var_minmax[index_right..., right] = minmax(var_minmax[index_right..., + right], var_left) + var_minmax[index_left..., left] = minmax(var_minmax[index_left..., left], + var_right) + end + end + + # Calc bounds at physical boundaries + for boundary in eachboundary(dg, cache) + element = cache.boundaries.neighbor_ids[boundary] + orientation = cache.boundaries.orientations[boundary] + neighbor_side = cache.boundaries.neighbor_sides[boundary] + + for i in eachnode(dg) + if neighbor_side == 2 # Element is on the right, boundary on the left + index = (1, i) + boundary_index = 1 + else # Element is on the left, boundary on the right + index = (nnodes(dg), i) + boundary_index = 2 + end + if orientation == 2 + index = reverse(index) + boundary_index += 2 + end + u_outer = get_boundary_outer_state(boundary_conditions[boundary_index], + cache, t, equations, dg, + index..., element) + var_outer = variable(u_outer, equations) + + var_minmax[index..., element] = minmax(var_minmax[index..., element], + var_outer) + end + end + + return nothing +end + ############################################################################### # Local minimum/maximum limiting @@ -232,6 +339,54 @@ end return nothing end +############################################################################## +# Local minimum limiting of specific entropy + +@inline function idp_spec_entropy!(alpha, limiter, u, t, dt, semi) + _, equations, dg, cache = mesh_equations_solver_cache(semi) + (; variable_bounds) = limiter.cache.subcell_limiter_coefficients + s_min = variable_bounds[:spec_entropy_min] + calc_bounds_onesided!(s_min, min, typemax, entropy_spec, u, t, semi) + + # Perform Newton's bisection method to find new alpha + @threaded for element in eachelement(dg, cache) + inverse_jacobian = cache.elements.inverse_jacobian[element] + for j in eachnode(dg), i in eachnode(dg) + u_local = get_node_vars(u, equations, dg, i, j, element) + newton_loops_alpha!(alpha, s_min[i, j, element], u_local, i, j, element, + entropy_spec, initial_check_entropy_spec, + final_check_nonnegative, inverse_jacobian, + dt, equations, dg, cache, limiter) + end + end + + return nothing +end + +############################################################################### +# Local maximum limiting of mathematical entropy + +@inline function idp_math_entropy!(alpha, limiter, u, t, dt, semi) + _, equations, dg, cache = mesh_equations_solver_cache(semi) + (; variable_bounds) = limiter.cache.subcell_limiter_coefficients + s_max = variable_bounds[:math_entropy_max] + calc_bounds_onesided!(s_max, max, typemin, entropy_math, u, t, semi) + + # Perform Newton's bisection method to find new alpha + @threaded for element in eachelement(dg, cache) + inverse_jacobian = cache.elements.inverse_jacobian[element] + for j in eachnode(dg), i in eachnode(dg) + u_local = get_node_vars(u, equations, dg, i, j, element) + newton_loops_alpha!(alpha, s_max[i, j, element], u_local, i, j, element, + entropy_math, initial_check_entropy_math, + final_check_nonnegative, inverse_jacobian, + dt, equations, dg, cache, limiter) + end + end + + return nothing +end + ############################################################################### # Global positivity limiting @@ -485,6 +640,14 @@ end end # Initial checks +@inline function initial_check_entropy_spec(bound, goal, newton_abstol) + goal <= max(newton_abstol, abs(bound) * newton_abstol) +end + +@inline function initial_check_entropy_math(bound, goal, newton_abstol) + goal >= -max(newton_abstol, abs(bound) * newton_abstol) +end + @inline initial_check_nonnegative(bound, goal, newton_abstol) = goal <= 0 # Goal and d(Goal)d(u) function @@ -493,7 +656,11 @@ end -dot(variable_derivative(variable, u, equations), dt * antidiffusive_flux) end -# Final checks +# Final check +@inline function final_check_standard(bound, goal, newton_abstol) + abs(goal) < max(newton_abstol, abs(bound) * newton_abstol) +end + @inline function final_check_nonnegative(bound, goal, newton_abstol) (goal <= eps()) && (goal > -max(newton_abstol, abs(bound) * newton_abstol)) end diff --git a/test/test_tree_2d_euler.jl b/test/test_tree_2d_euler.jl index 06b2da681cd..cb9be39b38d 100644 --- a/test/test_tree_2d_euler.jl +++ b/test/test_tree_2d_euler.jl @@ -317,16 +317,16 @@ end @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_blast_wave_sc_subcell_nonperiodic.jl"), l2=[ - 0.3517507570120483, - 0.19252291020146015, - 0.19249751956580294, - 0.618717827188004, + 0.3221133847370183, + 0.1798445067477554, + 0.17983199781537987, + 0.6136848341801259, ], linf=[ - 1.6699566795772216, - 1.3608007992899402, - 1.361864507190922, - 2.44022884092527, + 1.3433576692371516, + 1.1747340428151851, + 1.1744175397224497, + 2.4215287080587955, ], tspan=(0.0, 0.5), initial_refinement_level=4, @@ -371,16 +371,16 @@ end @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_sedov_blast_wave_sc_subcell.jl"), l2=[ - 0.4328635350273501, - 0.15011135840723572, - 0.15011135840723572, - 0.616129927549474, + 0.41446162783024115, + 0.14606475224088192, + 0.14607706259100364, + 0.6168364393501943, ], linf=[ - 1.6145297181778906, - 0.8614006163026988, - 0.8614006163026972, - 6.450225090647602, + 1.5717827990209163, + 0.7963182071213032, + 0.7956862000276138, + 6.456116489970984, ], tspan=(0.0, 1.0), initial_refinement_level=4, diff --git a/test/test_unit.jl b/test/test_unit.jl index 8b91f0875c4..16d94fabe9c 100644 --- a/test/test_unit.jl +++ b/test/test_unit.jl @@ -416,8 +416,8 @@ end indicator_hg = IndicatorHennemannGassner(1.0, 0.0, true, "variable", "cache") @test_nowarn show(stdout, indicator_hg) - limiter_idp = SubcellLimiterIDP(true, [1], true, [1], ["variable"], 0.1, "cache", 1, - (1.0, 1.0), 1.0) + limiter_idp = SubcellLimiterIDP(true, [1], true, [1], ["variable"], 0.1, true, true, + "cache", 1, (1.0, 1.0), 1.0) @test_nowarn show(stdout, limiter_idp) # TODO: TrixiShallowWater: move unit test From 927d783004ba319eb9e746c89d80087c3cab00a4 Mon Sep 17 00:00:00 2001 From: bennibolm Date: Fri, 22 Dec 2023 14:48:33 +0100 Subject: [PATCH 10/25] Update test errors after adding entropy limiting --- test/test_tree_2d_euler.jl | 32 ++++++++++++++++---------------- 1 file changed, 16 insertions(+), 16 deletions(-) diff --git a/test/test_tree_2d_euler.jl b/test/test_tree_2d_euler.jl index cb9be39b38d..1bb7cf89ea7 100644 --- a/test/test_tree_2d_euler.jl +++ b/test/test_tree_2d_euler.jl @@ -317,16 +317,16 @@ end @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_blast_wave_sc_subcell_nonperiodic.jl"), l2=[ - 0.3221133847370183, - 0.1798445067477554, - 0.17983199781537987, - 0.6136848341801259, + 0.32211638619244837, + 0.17984613698288662, + 0.1798341339489921, + 0.6136855864874968, ], linf=[ - 1.3433576692371516, - 1.1747340428151851, - 1.1744175397224497, - 2.4215287080587955, + 1.343528359596545, + 1.174708306029737, + 1.174481988232382, + 2.42152874599586, ], tspan=(0.0, 0.5), initial_refinement_level=4, @@ -371,16 +371,16 @@ end @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_sedov_blast_wave_sc_subcell.jl"), l2=[ - 0.41446162783024115, - 0.14606475224088192, - 0.14607706259100364, - 0.6168364393501943, + 0.4144617475432438, + 0.14606478267219714, + 0.1460770935963815, + 0.6168365992433479, ], linf=[ - 1.5717827990209163, - 0.7963182071213032, - 0.7956862000276138, - 6.456116489970984, + 1.5717829602926576, + 0.7963192532954594, + 0.7956878517936237, + 6.4561186367550825, ], tspan=(0.0, 1.0), initial_refinement_level=4, From e097c15fe8970fef3392bd6cb4a6269184e84524 Mon Sep 17 00:00:00 2001 From: bennibolm Date: Tue, 26 Dec 2023 16:58:01 +0100 Subject: [PATCH 11/25] Fix bug in entropy limiting --- src/solvers/dgsem_tree/subcell_limiters_2d.jl | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/src/solvers/dgsem_tree/subcell_limiters_2d.jl b/src/solvers/dgsem_tree/subcell_limiters_2d.jl index edc2a2182ad..ada4fd782f2 100644 --- a/src/solvers/dgsem_tree/subcell_limiters_2d.jl +++ b/src/solvers/dgsem_tree/subcell_limiters_2d.jl @@ -355,7 +355,7 @@ end u_local = get_node_vars(u, equations, dg, i, j, element) newton_loops_alpha!(alpha, s_min[i, j, element], u_local, i, j, element, entropy_spec, initial_check_entropy_spec, - final_check_nonnegative, inverse_jacobian, + final_check_standard, inverse_jacobian, dt, equations, dg, cache, limiter) end end @@ -379,7 +379,7 @@ end u_local = get_node_vars(u, equations, dg, i, j, element) newton_loops_alpha!(alpha, s_max[i, j, element], u_local, i, j, element, entropy_math, initial_check_entropy_math, - final_check_nonnegative, inverse_jacobian, + final_check_standard, inverse_jacobian, dt, equations, dg, cache, limiter) end end @@ -656,11 +656,13 @@ end -dot(variable_derivative(variable, u, equations), dt * antidiffusive_flux) end -# Final check +# Final checks +# final check for entropy limiting @inline function final_check_standard(bound, goal, newton_abstol) abs(goal) < max(newton_abstol, abs(bound) * newton_abstol) end +# final check for nonnegativity limiting @inline function final_check_nonnegative(bound, goal, newton_abstol) (goal <= eps()) && (goal > -max(newton_abstol, abs(bound) * newton_abstol)) end From 2e3ec712ea45f5c9a7891ac0561854b56c92cbdc Mon Sep 17 00:00:00 2001 From: bennibolm Date: Tue, 26 Dec 2023 17:29:12 +0100 Subject: [PATCH 12/25] Adapt estimated errors after bug fix --- test/test_tree_2d_euler.jl | 32 ++++++++++++++++---------------- 1 file changed, 16 insertions(+), 16 deletions(-) diff --git a/test/test_tree_2d_euler.jl b/test/test_tree_2d_euler.jl index 1bb7cf89ea7..c7e84cc7940 100644 --- a/test/test_tree_2d_euler.jl +++ b/test/test_tree_2d_euler.jl @@ -317,16 +317,16 @@ end @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_blast_wave_sc_subcell_nonperiodic.jl"), l2=[ - 0.32211638619244837, - 0.17984613698288662, - 0.1798341339489921, - 0.6136855864874968, + 0.32211935909459755, + 0.17984734129031393, + 0.17983517637561672, + 0.6136858672602447, ], linf=[ - 1.343528359596545, - 1.174708306029737, - 1.174481988232382, - 2.42152874599586, + 1.3435940880493509, + 1.1748248970276045, + 1.1745481442875036, + 2.421529617190895, ], tspan=(0.0, 0.5), initial_refinement_level=4, @@ -371,16 +371,16 @@ end @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_sedov_blast_wave_sc_subcell.jl"), l2=[ - 0.4144617475432438, - 0.14606478267219714, - 0.1460770935963815, - 0.6168365992433479, + 0.4144617128623457, + 0.14606474605745345, + 0.14607706120162497, + 0.616836505667339, ], linf=[ - 1.5717829602926576, - 0.7963192532954594, - 0.7956878517936237, - 6.4561186367550825, + 1.5717812331790788, + 0.7963178376801539, + 0.795685946979518, + 6.456116356402219, ], tspan=(0.0, 1.0), initial_refinement_level=4, From b33137cc41b592625f4b0be01a1aa57468063770 Mon Sep 17 00:00:00 2001 From: bennibolm Date: Wed, 31 Jan 2024 22:31:24 +0100 Subject: [PATCH 13/25] Remove doubled code --- src/callbacks_stage/subcell_bounds_check.jl | 5 ----- src/callbacks_stage/subcell_bounds_check_2d.jl | 6 ++++-- 2 files changed, 4 insertions(+), 7 deletions(-) diff --git a/src/callbacks_stage/subcell_bounds_check.jl b/src/callbacks_stage/subcell_bounds_check.jl index 4ca1de59a30..e3b8afb9391 100644 --- a/src/callbacks_stage/subcell_bounds_check.jl +++ b/src/callbacks_stage/subcell_bounds_check.jl @@ -164,11 +164,6 @@ end println(variable_string * ":\n- positivity: ", idp_bounds_delta_global[Symbol(variable_string, "_min")]) end - for variable in limiter.positivity_variables_nonlinear - variable_string = string(variable) - println(variable_string * ":\n- positivity: ", - idp_bounds_delta_global[Symbol(variable_string, "_min")]) - end end println("─"^100 * "\n") diff --git a/src/callbacks_stage/subcell_bounds_check_2d.jl b/src/callbacks_stage/subcell_bounds_check_2d.jl index 69d415b9799..61920565dde 100644 --- a/src/callbacks_stage/subcell_bounds_check_2d.jl +++ b/src/callbacks_stage/subcell_bounds_check_2d.jl @@ -49,7 +49,8 @@ @threaded for element in eachelement(solver, cache) deviation = deviation_threaded[stride_size * Threads.threadid()] for j in eachnode(solver), i in eachnode(solver) - s = entropy_spec(get_node_vars(u, equations, solver, i, j, element), equations) + s = entropy_spec(get_node_vars(u, equations, solver, i, j, element), + equations) deviation = max(deviation, variable_bounds[key][i, j, element] - s) end deviation_threaded[stride_size * Threads.threadid()] = deviation @@ -61,7 +62,8 @@ @threaded for element in eachelement(solver, cache) deviation = deviation_threaded[stride_size * Threads.threadid()] for j in eachnode(solver), i in eachnode(solver) - s = entropy_math(get_node_vars(u, equations, solver, i, j, element), equations) + s = entropy_math(get_node_vars(u, equations, solver, i, j, element), + equations) deviation = max(deviation, s - variable_bounds[key][i, j, element]) end deviation_threaded[stride_size * Threads.threadid()] = deviation From 9cf4b8e7045204ee7f57c25054ce94b4b27f531c Mon Sep 17 00:00:00 2001 From: bennibolm Date: Wed, 31 Jan 2024 22:34:54 +0100 Subject: [PATCH 14/25] Rename function --- src/equations/compressible_euler_2d.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/equations/compressible_euler_2d.jl b/src/equations/compressible_euler_2d.jl index 1ec9576dd19..2cb7fb3d028 100644 --- a/src/equations/compressible_euler_2d.jl +++ b/src/equations/compressible_euler_2d.jl @@ -1721,8 +1721,8 @@ end end # Transformation from conservative variables u to d(s)/d(u) -@inline function variable_derivative(::typeof(entropy_math), - u, equations::CompressibleEulerEquations2D) +@inline function gradient_conservative(::typeof(entropy_math), + u, equations::CompressibleEulerEquations2D) return cons2entropy(u, equations) end @@ -1740,8 +1740,8 @@ end end # Transformation from conservative variables u to d(s)/d(u) -@inline function variable_derivative(::typeof(entropy_spec), - u, equations::CompressibleEulerEquations2D) +@inline function gradient_conservative(::typeof(entropy_spec), + u, equations::CompressibleEulerEquations2D) return cons2entropy_spec(u, equations) end From 37b81172b25df9a8b56f040f1020d2fa61cdd6f2 Mon Sep 17 00:00:00 2001 From: Benjamin Bolm <74359358+bennibolm@users.noreply.github.com> Date: Fri, 1 Mar 2024 13:58:36 +0100 Subject: [PATCH 15/25] Generalize one-sided limiting (#125) * Start to align both entropy limiters * Adapt calc_bounds_onesided! * Add wrapper function for entropy limiting * Rename keys in Dict * use variable and bound_function as parameter * Use same function for both entropies * First working version of general onesided limiting * Rename minmax limiting to twosided limiting * Update comment * Clean up default vector * Last stuff * Fix unit test * fmt * Fix tests * Correct order * Rework docstring * Rename operator to min_or_max * Call initial check with min_or_max * fmt * Implement suggestions * Remove type stuff * Fix allocations due to non-specialized routine --- ...euler_blast_wave_sc_subcell_nonperiodic.jl | 5 +- ...lixir_euler_sedov_blast_wave_sc_subcell.jl | 5 +- ...ck_bubble_shockcapturing_subcell_minmax.jl | 4 +- src/callbacks_stage/subcell_bounds_check.jl | 41 ++--- .../subcell_bounds_check_2d.jl | 64 ++++--- src/solvers/dgsem_tree/subcell_limiters.jl | 119 +++++++------ src/solvers/dgsem_tree/subcell_limiters_2d.jl | 163 ++++++++---------- test/test_unit.jl | 5 +- 8 files changed, 204 insertions(+), 202 deletions(-) diff --git a/examples/tree_2d_dgsem/elixir_euler_blast_wave_sc_subcell_nonperiodic.jl b/examples/tree_2d_dgsem/elixir_euler_blast_wave_sc_subcell_nonperiodic.jl index 6ff381b896c..00d3c69f2e6 100644 --- a/examples/tree_2d_dgsem/elixir_euler_blast_wave_sc_subcell_nonperiodic.jl +++ b/examples/tree_2d_dgsem/elixir_euler_blast_wave_sc_subcell_nonperiodic.jl @@ -41,8 +41,9 @@ surface_flux = flux_lax_friedrichs volume_flux = flux_ranocha basis = LobattoLegendreBasis(3) limiter_idp = SubcellLimiterIDP(equations, basis; - local_minmax_variables_cons = ["rho"], - math_entropy = true) + local_twosided_variables_cons = ["rho"], + local_onesided_variables_nonlinear = [(Trixi.entropy_math, + max)]) volume_integral = VolumeIntegralSubcellLimiting(limiter_idp; volume_flux_dg = volume_flux, volume_flux_fv = surface_flux) diff --git a/examples/tree_2d_dgsem/elixir_euler_sedov_blast_wave_sc_subcell.jl b/examples/tree_2d_dgsem/elixir_euler_sedov_blast_wave_sc_subcell.jl index 64cae963158..d3bb6a819b5 100644 --- a/examples/tree_2d_dgsem/elixir_euler_sedov_blast_wave_sc_subcell.jl +++ b/examples/tree_2d_dgsem/elixir_euler_sedov_blast_wave_sc_subcell.jl @@ -42,8 +42,9 @@ surface_flux = flux_lax_friedrichs volume_flux = flux_chandrashekar basis = LobattoLegendreBasis(3) limiter_idp = SubcellLimiterIDP(equations, basis; - local_minmax_variables_cons = ["rho"], - spec_entropy = true) + local_twosided_variables_cons = ["rho"], + local_onesided_variables_nonlinear = [(Trixi.entropy_spec, + min)]) volume_integral = VolumeIntegralSubcellLimiting(limiter_idp; volume_flux_dg = volume_flux, volume_flux_fv = surface_flux) diff --git a/examples/tree_2d_dgsem/elixir_eulermulti_shock_bubble_shockcapturing_subcell_minmax.jl b/examples/tree_2d_dgsem/elixir_eulermulti_shock_bubble_shockcapturing_subcell_minmax.jl index 4b606502ebe..b2d49ecbd48 100644 --- a/examples/tree_2d_dgsem/elixir_eulermulti_shock_bubble_shockcapturing_subcell_minmax.jl +++ b/examples/tree_2d_dgsem/elixir_eulermulti_shock_bubble_shockcapturing_subcell_minmax.jl @@ -86,8 +86,8 @@ volume_flux = flux_ranocha basis = LobattoLegendreBasis(3) limiter_idp = SubcellLimiterIDP(equations, basis; - local_minmax_variables_cons = ["rho" * string(i) - for i in eachcomponent(equations)]) + local_twosided_variables_cons = ["rho" * string(i) + for i in eachcomponent(equations)]) volume_integral = VolumeIntegralSubcellLimiting(limiter_idp; volume_flux_dg = volume_flux, volume_flux_fv = surface_flux) diff --git a/src/callbacks_stage/subcell_bounds_check.jl b/src/callbacks_stage/subcell_bounds_check.jl index e3b8afb9391..ba193ab2997 100644 --- a/src/callbacks_stage/subcell_bounds_check.jl +++ b/src/callbacks_stage/subcell_bounds_check.jl @@ -77,28 +77,27 @@ function init_callback(callback::BoundsCheckCallback, semi, limiter::SubcellLimi return nothing end - (; local_minmax, positivity, spec_entropy, math_entropy) = limiter + (; local_twosided, positivity, local_onesided) = 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 local_minmax - for v in limiter.local_minmax_variables_cons + if local_twosided + for v in limiter.local_twosided_variables_cons variable_string = string(variables[v]) print(f, ", " * variable_string * "_min, " * variable_string * "_max") end end - if spec_entropy - print(f, ", specEntr_min") - end - if math_entropy - print(f, ", mathEntr_max") + if local_onesided + for (variable, min_or_max) in limiter.local_onesided_variables_nonlinear + print(f, ", " * string(variable) * "_" * string(min_or_max)) + end end if positivity for v in limiter.positivity_variables_cons - if v in limiter.local_minmax_variables_cons + if v in limiter.local_twosided_variables_cons continue end print(f, ", " * string(variables[v]) * "_min") @@ -126,15 +125,15 @@ end @inline function finalize_callback(callback::BoundsCheckCallback, semi, limiter::SubcellLimiterIDP) - (; local_minmax, positivity, spec_entropy, math_entropy) = limiter + (; local_twosided, positivity, local_onesided) = limiter (; idp_bounds_delta_global) = limiter.cache variables = varnames(cons2cons, semi.equations) println("─"^100) println("Maximum deviation from bounds:") println("─"^100) - if local_minmax - for v in limiter.local_minmax_variables_cons + if local_twosided + for v in limiter.local_twosided_variables_cons v_string = string(v) println("$(variables[v]):") println("- lower bound: ", @@ -143,17 +142,19 @@ end idp_bounds_delta_global[Symbol(v_string, "_max")]) end end - if spec_entropy - println("spec. entropy:\n- lower bound: ", - idp_bounds_delta_global[:spec_entropy_min]) - end - if math_entropy - println("math. entropy:\n- upper bound: ", - idp_bounds_delta_global[:math_entropy_max]) + if local_onesided + for (variable, min_or_max) in limiter.local_onesided_variables_nonlinear + variable_string = string(variable) + minmax_string = string(min_or_max) + println("$variable_string:") + println("- $minmax_string bound: ", + idp_bounds_delta_global[Symbol(variable_string, "_", + minmax_string)]) + end end if positivity for v in limiter.positivity_variables_cons - if v in limiter.local_minmax_variables_cons + if v in limiter.local_twosided_variables_cons continue end println(string(variables[v]) * ":\n- positivity: ", diff --git a/src/callbacks_stage/subcell_bounds_check_2d.jl b/src/callbacks_stage/subcell_bounds_check_2d.jl index 61920565dde..4b8fb34c6b5 100644 --- a/src/callbacks_stage/subcell_bounds_check_2d.jl +++ b/src/callbacks_stage/subcell_bounds_check_2d.jl @@ -8,7 +8,7 @@ @inline function check_bounds(u, mesh::AbstractMesh{2}, equations, solver, cache, limiter::SubcellLimiterIDP, time, iter, output_directory, save_errors) - (; local_minmax, positivity, spec_entropy, math_entropy) = solver.volume_integral.limiter + (; local_twosided, positivity, local_onesided) = solver.volume_integral.limiter (; variable_bounds) = limiter.cache.subcell_limiter_coefficients (; idp_bounds_delta_local, idp_bounds_delta_global) = limiter.cache @@ -21,8 +21,8 @@ # Since there are no processors with caches over 128B, we use `n = 128B / size(uEltype)` stride_size = div(128, sizeof(eltype(u))) # = n - if local_minmax - for v in limiter.local_minmax_variables_cons + if local_twosided + for v in limiter.local_twosided_variables_cons v_string = string(v) key_min = Symbol(v_string, "_min") key_max = Symbol(v_string, "_max") @@ -43,35 +43,30 @@ end end end - if spec_entropy - key = :spec_entropy_min - deviation_threaded = idp_bounds_delta_local[key] - @threaded for element in eachelement(solver, cache) - deviation = deviation_threaded[stride_size * Threads.threadid()] - for j in eachnode(solver), i in eachnode(solver) - s = entropy_spec(get_node_vars(u, equations, solver, i, j, element), - equations) - deviation = max(deviation, variable_bounds[key][i, j, element] - s) - end - deviation_threaded[stride_size * Threads.threadid()] = deviation - end - end - if math_entropy - key = :math_entropy_max - deviation_threaded = idp_bounds_delta_local[key] - @threaded for element in eachelement(solver, cache) - deviation = deviation_threaded[stride_size * Threads.threadid()] - for j in eachnode(solver), i in eachnode(solver) - s = entropy_math(get_node_vars(u, equations, solver, i, j, element), + if local_onesided + foreach(limiter.local_onesided_variables_nonlinear) do (variable, min_or_max) + key = Symbol(string(variable), "_", string(min_or_max)) + deviation_threaded = idp_bounds_delta_local[key] + @threaded for element in eachelement(solver, cache) + deviation = deviation_threaded[stride_size * Threads.threadid()] + for j in eachnode(solver), i in eachnode(solver) + v = variable(get_node_vars(u, equations, solver, i, j, element), equations) - deviation = max(deviation, s - variable_bounds[key][i, j, element]) + if min_or_max === max + deviation = max(deviation, + v - variable_bounds[key][i, j, element]) + else # min_or_max === min + deviation = max(deviation, + variable_bounds[key][i, j, element] - v) + end + end + deviation_threaded[stride_size * Threads.threadid()] = deviation end - deviation_threaded[stride_size * Threads.threadid()] = deviation end end if positivity for v in limiter.positivity_variables_cons - if v in limiter.local_minmax_variables_cons + if v in limiter.local_twosided_variables_cons continue end key = Symbol(string(v), "_min") @@ -115,8 +110,8 @@ # Print to output file open("$output_directory/deviations.txt", "a") do f print(f, iter, ", ", time) - if local_minmax - for v in limiter.local_minmax_variables_cons + if local_twosided + for v in limiter.local_twosided_variables_cons v_string = string(v) print(f, ", ", idp_bounds_delta_local[Symbol(v_string, "_min")][stride_size], @@ -124,15 +119,16 @@ idp_bounds_delta_local[Symbol(v_string, "_max")][stride_size]) end end - if spec_entropy - print(f, ", ", idp_bounds_delta_local[:spec_entropy_min][stride_size]) - end - if math_entropy - print(f, ", ", idp_bounds_delta_local[:math_entropy_max][stride_size]) + if local_onesided + for (variable, min_or_max) in limiter.local_onesided_variables_nonlinear + print(f, ", ", + idp_bounds_delta_local[Symbol(string(variable), "_", + string(min_or_max))][stride_size]) + end end if positivity for v in limiter.positivity_variables_cons - if v in limiter.local_minmax_variables_cons + if v in limiter.local_twosided_variables_cons continue end print(f, ", ", diff --git a/src/solvers/dgsem_tree/subcell_limiters.jl b/src/solvers/dgsem_tree/subcell_limiters.jl index bd865b8e1e4..7b06149148d 100644 --- a/src/solvers/dgsem_tree/subcell_limiters.jl +++ b/src/solvers/dgsem_tree/subcell_limiters.jl @@ -14,30 +14,36 @@ end """ SubcellLimiterIDP(equations::AbstractEquations, basis; - local_minmax_variables_cons = String[], + local_twosided_variables_cons = String[], positivity_variables_cons = String[], positivity_variables_nonlinear = [], positivity_correction_factor = 0.1, - spec_entropy = false, - math_entropy = false, + local_onesided_variables_nonlinear = [], max_iterations_newton = 10, newton_tolerances = (1.0e-12, 1.0e-14), gamma_constant_newton = 2 * ndims(equations)) Subcell invariant domain preserving (IDP) limiting used with [`VolumeIntegralSubcellLimiting`](@ref) including: -- Local maximum/minimum Zalesak-type limiting for conservative variables (`local_minmax_variables_cons`) +- Local two-sided Zalesak-type limiting for conservative variables (`local_twosided_variables_cons`) - Positivity limiting for conservative variables (`positivity_variables_cons`) and nonlinear variables (`positivity_variables_nonlinear`) -- One-sided limiting for specific and mathematical entropy (`spec_entropy`, `math_entropy`) +- Local one-sided limiting for nonlinear variables, e.g. `entropy_spec` and `entropy_math` +with `local_onesided_variables_nonlinear` -Conservative variables to be limited are passed as a vector of strings, e.g. `local_minmax_variables_cons = ["rho"]` -and `positivity_variables_cons = ["rho"]`. For nonlinear variables the specific functions are -passed in a vector, e.g. `positivity_variables_nonlinear = [pressure]`. +To use these three limiting options use the following structure: + +***Conservative variables*** to be limited are passed as a vector of strings, e.g. +`local_twosided_variables_cons = ["rho"]` and `positivity_variables_cons = ["rho"]`. +For ***nonlinear variables***, the wanted variable functions are passed within a vector: To ensure +positivity use a plain vector including the desired variables, e.g. `positivity_variables_nonlinear = [pressure]`. +For local one-sided limiting pass the variable function combined with the requested bound +(`min` or `max`) as a tuple. For instance, to impose a lower local bound on the specific entropy +use `local_onesided_variables_nonlinear = [(Trixi.entropy_spec, min)]`. The bounds are calculated using the low-order FV solution. The positivity limiter uses `positivity_correction_factor` such that `u^new >= positivity_correction_factor * u^FV`. -The limiting of nonlinear variables uses a Newton-bisection method with a maximum of +Local and global limiting of nonlinear variables uses a Newton-bisection method with a maximum of `max_iterations_newton` iterations, relative and absolute tolerances of `newton_tolerances` and a provisional update constant `gamma_constant_newton` (`gamma_constant_newton>=2*d`, where `d = #dimensions`). See equation (20) of Pazner (2020) and equation (30) of Rueda-Ramírez et al. (2022). @@ -58,16 +64,17 @@ where `d = #dimensions`). See equation (20) of Pazner (2020) and equation (30) o !!! warning "Experimental implementation" This is an experimental feature and may change in future releases. """ -struct SubcellLimiterIDP{RealT <: Real, LimitingVariablesNonlinear, Cache} <: +struct SubcellLimiterIDP{RealT <: Real, LimitingVariablesNonlinear, + LimitingOnesidedVariablesNonlinear, Cache} <: AbstractSubcellLimiter - local_minmax::Bool - local_minmax_variables_cons::Vector{Int} # Local mininum/maximum principles for conservative variables + local_twosided::Bool + local_twosided_variables_cons::Vector{Int} # Local two-sided limiting for conservative variables positivity::Bool positivity_variables_cons::Vector{Int} # Positivity for conservative variables positivity_variables_nonlinear::LimitingVariablesNonlinear # Positivity for nonlinear variables positivity_correction_factor::RealT - spec_entropy::Bool - math_entropy::Bool + local_onesided::Bool + local_onesided_variables_nonlinear::LimitingOnesidedVariablesNonlinear # Local one-sided limiting for nonlinear variables cache::Cache max_iterations_newton::Int newton_tolerances::Tuple{RealT, RealT} # Relative and absolute tolerances for Newton's method @@ -76,43 +83,56 @@ end # this method is used when the limiter is constructed as for shock-capturing volume integrals function SubcellLimiterIDP(equations::AbstractEquations, basis; - local_minmax_variables_cons = String[], + local_twosided_variables_cons = String[], positivity_variables_cons = String[], positivity_variables_nonlinear = [], positivity_correction_factor = 0.1, - spec_entropy = false, - math_entropy = false, + local_onesided_variables_nonlinear = [], max_iterations_newton = 10, newton_tolerances = (1.0e-12, 1.0e-14), gamma_constant_newton = 2 * ndims(equations)) - local_minmax = (length(local_minmax_variables_cons) > 0) + local_twosided = (length(local_twosided_variables_cons) > 0) + local_onesided = (length(local_onesided_variables_nonlinear) > 0) positivity = (length(positivity_variables_cons) + length(positivity_variables_nonlinear) > 0) - if math_entropy && spec_entropy - error("Only one of the two can be selected: math_entropy/spec_entropy") + + # When passing `min` or `max` in the elixir, the specific function of Base is used. + # To speed up the simulation, we replace it with `Trixi.min` and `Trixi.max` respectively. + local_onesided_variables_nonlinear_ = Tuple{Function, Function}[] + for (variable, min_or_max) in local_onesided_variables_nonlinear + if min_or_max === Base.max + push!(local_onesided_variables_nonlinear_, tuple(variable, max)) + elseif min_or_max === Base.min + push!(local_onesided_variables_nonlinear_, tuple(variable, min)) + elseif min_or_max === Trixi.max || min_or_max === Trixi.min + push!(local_onesided_variables_nonlinear_, tuple(variable, min_or_max)) + else + error("Parameter $min_or_max is not a valid input. Use `max` or `min` instead.") + end end + local_onesided_variables_nonlinear_ = Tuple(local_onesided_variables_nonlinear_) - local_minmax_variables_cons_ = get_variable_index.(local_minmax_variables_cons, - equations) + local_twosided_variables_cons_ = get_variable_index.(local_twosided_variables_cons, + equations) positivity_variables_cons_ = get_variable_index.(positivity_variables_cons, equations) bound_keys = () - if local_minmax - for v in local_minmax_variables_cons_ + if local_twosided + for v in local_twosided_variables_cons_ v_string = string(v) bound_keys = (bound_keys..., Symbol(v_string, "_min"), Symbol(v_string, "_max")) end end - if spec_entropy - bound_keys = (bound_keys..., :spec_entropy_min) - end - if math_entropy - bound_keys = (bound_keys..., :math_entropy_max) + if local_onesided + for (variable, min_or_max) in local_onesided_variables_nonlinear_ + bound_keys = (bound_keys..., + Symbol(string(variable), "_", string(min_or_max))) + end end for v in positivity_variables_cons_ - if !(v in local_minmax_variables_cons_) + if !(v in local_twosided_variables_cons_) bound_keys = (bound_keys..., Symbol(string(v), "_min")) end end @@ -124,11 +144,13 @@ function SubcellLimiterIDP(equations::AbstractEquations, basis; SubcellLimiterIDP{typeof(positivity_correction_factor), typeof(positivity_variables_nonlinear), - typeof(cache)}(local_minmax, local_minmax_variables_cons_, + typeof(local_onesided_variables_nonlinear_), + typeof(cache)}(local_twosided, local_twosided_variables_cons_, positivity, positivity_variables_cons_, positivity_variables_nonlinear, positivity_correction_factor, - spec_entropy, math_entropy, + local_onesided, + local_onesided_variables_nonlinear_, cache, max_iterations_newton, newton_tolerances, gamma_constant_newton) @@ -136,24 +158,21 @@ end function Base.show(io::IO, limiter::SubcellLimiterIDP) @nospecialize limiter # reduce precompilation time - (; local_minmax, positivity, spec_entropy, math_entropy) = limiter + (; local_twosided, positivity, local_onesided) = limiter print(io, "SubcellLimiterIDP(") - if !(local_minmax || positivity || spec_entropy || math_entropy) + if !(local_twosided || positivity || local_onesided) print(io, "No limiter selected => pure DG method") else features = String[] - if local_minmax + if local_twosided push!(features, "local min/max") end if positivity push!(features, "positivity") end - if spec_entropy - push!(features, "specific entropy") - end - if math_entropy - push!(features, "mathematical entropy") + if local_onesided + push!(features, "local onesided") end join(io, features, ", ") print(io, "Limiter=($features), ") @@ -164,19 +183,19 @@ end function Base.show(io::IO, ::MIME"text/plain", limiter::SubcellLimiterIDP) @nospecialize limiter # reduce precompilation time - (; local_minmax, positivity, spec_entropy, math_entropy) = limiter + (; local_twosided, positivity, local_onesided) = limiter if get(io, :compact, false) show(io, limiter) else - if !(local_minmax || positivity || spec_entropy || math_entropy) + if !(local_twosided || positivity || local_onesided) setup = ["Limiter" => "No limiter selected => pure DG method"] else setup = ["Limiter" => ""] - if local_minmax + if local_twosided setup = [ setup..., - "" => "Local maximum/minimum limiting for conservative variables $(limiter.local_minmax_variables_cons)", + "" => "Local two-sided limiting for conservative variables $(limiter.local_twosided_variables_cons)", ] end if positivity @@ -187,14 +206,10 @@ function Base.show(io::IO, ::MIME"text/plain", limiter::SubcellLimiterIDP) "" => "- with positivity correction factor = $(limiter.positivity_correction_factor)", ] end - if spec_entropy - setup = [setup..., "" => "Local minimum limiting for specific entropy"] - end - if math_entropy - setup = [ - setup..., - "" => "Local maximum limiting for mathematical entropy", - ] + if local_onesided + for (variable, min_or_max) in limiter.local_onesided_variables_nonlinear + setup = [setup..., "" => "Local $min_or_max limiting for $variable"] + end end setup = [ setup..., diff --git a/src/solvers/dgsem_tree/subcell_limiters_2d.jl b/src/solvers/dgsem_tree/subcell_limiters_2d.jl index 3d831e4fa39..e820bed667c 100644 --- a/src/solvers/dgsem_tree/subcell_limiters_2d.jl +++ b/src/solvers/dgsem_tree/subcell_limiters_2d.jl @@ -44,20 +44,16 @@ function (limiter::SubcellLimiterIDP)(u::AbstractArray{<:Any, 4}, semi, dg::DGSE # TODO: Do not abuse `reset_du!` but maybe implement a generic `set_zero!` @trixi_timeit timer() "reset alpha" reset_du!(alpha, dg, semi.cache) - if limiter.local_minmax - @trixi_timeit timer() "local min/max limiting" idp_local_minmax!(alpha, limiter, - u, t, dt, semi) + if limiter.local_twosided + @trixi_timeit timer() "local twosided" idp_local_twosided!(alpha, limiter, + u, t, dt, semi) end if limiter.positivity @trixi_timeit timer() "positivity" idp_positivity!(alpha, limiter, u, dt, semi) end - if limiter.spec_entropy - @trixi_timeit timer() "spec_entropy" idp_spec_entropy!(alpha, limiter, u, t, dt, - semi) - end - if limiter.math_entropy - @trixi_timeit timer() "math_entropy" idp_math_entropy!(alpha, limiter, u, t, dt, - semi) + if limiter.local_onesided + @trixi_timeit timer() "local onesided" idp_local_onesided!(alpha, limiter, + u, t, dt, semi) end # Calculate alpha1 and alpha2 @@ -179,44 +175,48 @@ end return nothing end -@inline function calc_bounds_onesided!(var_minmax, minmax, typeminmax, variable, u, t, - semi) +@inline function calc_bounds_onesided!(var_minmax, min_or_max, variable, u, t, semi) mesh, equations, dg, cache = mesh_equations_solver_cache(semi) # Calc bounds inside elements @threaded for element in eachelement(dg, cache) + # Reset bounds for j in eachnode(dg), i in eachnode(dg) - var_minmax[i, j, element] = typeminmax(eltype(var_minmax)) + if min_or_max === max + var_minmax[i, j, element] = typemin(eltype(var_minmax)) + else + var_minmax[i, j, element] = typemax(eltype(var_minmax)) + end end # Calculate bounds at Gauss-Lobatto nodes using u for j in eachnode(dg), i in eachnode(dg) var = variable(get_node_vars(u, equations, dg, i, j, element), equations) - var_minmax[i, j, element] = minmax(var_minmax[i, j, element], var) + var_minmax[i, j, element] = min_or_max(var_minmax[i, j, element], var) if i > 1 - var_minmax[i - 1, j, element] = minmax(var_minmax[i - 1, j, element], - var) + var_minmax[i - 1, j, element] = min_or_max(var_minmax[i - 1, j, + element], var) end if i < nnodes(dg) - var_minmax[i + 1, j, element] = minmax(var_minmax[i + 1, j, element], - var) + var_minmax[i + 1, j, element] = min_or_max(var_minmax[i + 1, j, + element], var) end if j > 1 - var_minmax[i, j - 1, element] = minmax(var_minmax[i, j - 1, element], - var) + var_minmax[i, j - 1, element] = min_or_max(var_minmax[i, j - 1, + element], var) end if j < nnodes(dg) - var_minmax[i, j + 1, element] = minmax(var_minmax[i, j + 1, element], - var) + var_minmax[i, j + 1, element] = min_or_max(var_minmax[i, j + 1, + element], var) end end end # Values at element boundary - calc_bounds_onesided_interface!(var_minmax, minmax, variable, u, t, semi, mesh) + calc_bounds_onesided_interface!(var_minmax, min_or_max, variable, u, t, semi, mesh) end -@inline function calc_bounds_onesided_interface!(var_minmax, minmax, variable, u, t, +@inline function calc_bounds_onesided_interface!(var_minmax, min_or_max, variable, u, t, semi, mesh::TreeMesh2D) _, equations, dg, cache = mesh_equations_solver_cache(semi) (; boundary_conditions) = semi @@ -240,10 +240,10 @@ end var_right = variable(get_node_vars(u, equations, dg, index_right..., right), equations) - var_minmax[index_right..., right] = minmax(var_minmax[index_right..., - right], var_left) - var_minmax[index_left..., left] = minmax(var_minmax[index_left..., left], - var_right) + var_minmax[index_right..., right] = min_or_max(var_minmax[index_right..., + right], var_left) + var_minmax[index_left..., left] = min_or_max(var_minmax[index_left..., + left], var_right) end end @@ -270,8 +270,8 @@ end index..., element) var_outer = variable(u_outer, equations) - var_minmax[index..., element] = minmax(var_minmax[index..., element], - var_outer) + var_minmax[index..., element] = min_or_max(var_minmax[index..., element], + var_outer) end end @@ -279,17 +279,17 @@ end end ############################################################################### -# Local minimum/maximum limiting +# Local two-sided limiting of conservative variables -@inline function idp_local_minmax!(alpha, limiter, u, t, dt, semi) - for variable in limiter.local_minmax_variables_cons - idp_local_minmax!(alpha, limiter, u, t, dt, semi, variable) +@inline function idp_local_twosided!(alpha, limiter, u, t, dt, semi) + for variable in limiter.local_twosided_variables_cons + idp_local_twosided!(alpha, limiter, u, t, dt, semi, variable) end return nothing end -@inline function idp_local_minmax!(alpha, limiter, u, t, dt, semi, variable) +@inline function idp_local_twosided!(alpha, limiter, u, t, dt, semi, variable) _, _, dg, cache = mesh_equations_solver_cache(semi) (; antidiffusive_flux1_L, antidiffusive_flux2_L, antidiffusive_flux1_R, antidiffusive_flux2_R) = cache.antidiffusive_fluxes (; inverse_weights) = dg.basis @@ -351,46 +351,32 @@ end end ############################################################################## -# Local minimum limiting of specific entropy - -@inline function idp_spec_entropy!(alpha, limiter, u, t, dt, semi) - _, equations, dg, cache = mesh_equations_solver_cache(semi) - (; variable_bounds) = limiter.cache.subcell_limiter_coefficients - s_min = variable_bounds[:spec_entropy_min] - calc_bounds_onesided!(s_min, min, typemax, entropy_spec, u, t, semi) +# Local one-sided limiting of nonlinear variables - # Perform Newton's bisection method to find new alpha - @threaded for element in eachelement(dg, cache) - inverse_jacobian = cache.elements.inverse_jacobian[element] - for j in eachnode(dg), i in eachnode(dg) - u_local = get_node_vars(u, equations, dg, i, j, element) - newton_loops_alpha!(alpha, s_min[i, j, element], u_local, i, j, element, - entropy_spec, initial_check_entropy_spec_newton_idp, - final_check_standard_newton_idp, inverse_jacobian, - dt, equations, dg, cache, limiter) - end +@inline function idp_local_onesided!(alpha, limiter, u, t, dt, semi) + foreach(limiter.local_onesided_variables_nonlinear) do (variable, min_or_max) + idp_local_onesided!(alpha, limiter, u, t, dt, semi, variable, min_or_max) end return nothing end -############################################################################### -# Local maximum limiting of mathematical entropy - -@inline function idp_math_entropy!(alpha, limiter, u, t, dt, semi) +@inline function idp_local_onesided!(alpha, limiter, u, t, dt, semi, variable, + min_or_max) _, equations, dg, cache = mesh_equations_solver_cache(semi) (; variable_bounds) = limiter.cache.subcell_limiter_coefficients - s_max = variable_bounds[:math_entropy_max] - calc_bounds_onesided!(s_max, max, typemin, entropy_math, u, t, semi) + var_minmax = variable_bounds[Symbol(string(variable), "_", string(min_or_max))] + calc_bounds_onesided!(var_minmax, min_or_max, variable, u, t, semi) # Perform Newton's bisection method to find new alpha @threaded for element in eachelement(dg, cache) inverse_jacobian = cache.elements.inverse_jacobian[element] for j in eachnode(dg), i in eachnode(dg) u_local = get_node_vars(u, equations, dg, i, j, element) - newton_loops_alpha!(alpha, s_max[i, j, element], u_local, i, j, element, - entropy_math, initial_check_entropy_math_newton_idp, - final_check_standard_newton_idp, inverse_jacobian, + newton_loops_alpha!(alpha, var_minmax[i, j, element], u_local, + i, j, element, variable, min_or_max, + initial_check_local_onesided_newton_idp, + final_check_local_onesided_newton_idp, inverse_jacobian, dt, equations, dg, cache, limiter) end end @@ -445,8 +431,8 @@ end end # Compute bound - if limiter.local_minmax && - variable in limiter.local_minmax_variables_cons && + if limiter.local_twosided && + variable in limiter.local_twosided_variables_cons && var_min[i, j, element] >= positivity_correction_factor * var # Local limiting is more restrictive that positivity limiting # => Skip positivity limiting for this node @@ -508,7 +494,7 @@ end # Perform Newton's bisection method to find new alpha newton_loops_alpha!(alpha, var_min[i, j, element], u_local, i, j, element, - variable, initial_check_nonnegative_newton_idp, + variable, min, initial_check_nonnegative_newton_idp, final_check_nonnegative_newton_idp, inverse_jacobian, dt, equations, dg, cache, limiter) end @@ -518,8 +504,9 @@ end end @inline function newton_loops_alpha!(alpha, bound, u, i, j, element, variable, - initial_check, final_check, inverse_jacobian, dt, - equations, dg, cache, limiter) + min_or_max, initial_check, final_check, + inverse_jacobian, dt, equations, dg, cache, + limiter) (; inverse_weights) = dg.basis (; antidiffusive_flux1_L, antidiffusive_flux2_L, antidiffusive_flux1_R, antidiffusive_flux2_R) = cache.antidiffusive_fluxes @@ -529,37 +516,38 @@ end antidiffusive_flux = gamma_constant_newton * inverse_jacobian * inverse_weights[i] * get_node_vars(antidiffusive_flux1_R, equations, dg, i, j, element) - newton_loop!(alpha, bound, u, i, j, element, variable, initial_check, final_check, - equations, dt, limiter, antidiffusive_flux) + newton_loop!(alpha, bound, u, i, j, element, variable, min_or_max, initial_check, + final_check, equations, dt, limiter, antidiffusive_flux) # positive xi direction antidiffusive_flux = -gamma_constant_newton * inverse_jacobian * inverse_weights[i] * get_node_vars(antidiffusive_flux1_L, equations, dg, i + 1, j, element) - newton_loop!(alpha, bound, u, i, j, element, variable, initial_check, final_check, - equations, dt, limiter, antidiffusive_flux) + newton_loop!(alpha, bound, u, i, j, element, variable, min_or_max, initial_check, + final_check, equations, dt, limiter, antidiffusive_flux) # negative eta direction antidiffusive_flux = gamma_constant_newton * inverse_jacobian * inverse_weights[j] * get_node_vars(antidiffusive_flux2_R, equations, dg, i, j, element) - newton_loop!(alpha, bound, u, i, j, element, variable, initial_check, final_check, - equations, dt, limiter, antidiffusive_flux) + newton_loop!(alpha, bound, u, i, j, element, variable, min_or_max, initial_check, + final_check, equations, dt, limiter, antidiffusive_flux) # positive eta direction antidiffusive_flux = -gamma_constant_newton * inverse_jacobian * inverse_weights[j] * get_node_vars(antidiffusive_flux2_L, equations, dg, i, j + 1, element) - newton_loop!(alpha, bound, u, i, j, element, variable, initial_check, final_check, - equations, dt, limiter, antidiffusive_flux) + newton_loop!(alpha, bound, u, i, j, element, variable, min_or_max, initial_check, + final_check, equations, dt, limiter, antidiffusive_flux) return nothing end -@inline function newton_loop!(alpha, bound, u, i, j, element, variable, initial_check, - final_check, equations, dt, limiter, antidiffusive_flux) +@inline function newton_loop!(alpha, bound, u, i, j, element, variable, min_or_max, + initial_check, final_check, equations, dt, limiter, + antidiffusive_flux) newton_reltol, newton_abstol = limiter.newton_tolerances beta = 1 - alpha[i, j, element] @@ -573,7 +561,7 @@ end if isvalid(u_curr, equations) goal = goal_function_newton_idp(variable, bound, u_curr, equations) - initial_check(bound, goal, newton_abstol) && return nothing + initial_check(min_or_max, bound, goal, newton_abstol) && return nothing end # Newton iterations @@ -608,7 +596,7 @@ end # Check new beta for condition and update bounds goal = goal_function_newton_idp(variable, bound, u_curr, equations) - if initial_check(bound, goal, newton_abstol) + if initial_check(min_or_max, bound, goal, newton_abstol) # New beta fulfills condition beta_L = beta else @@ -641,26 +629,25 @@ end end new_alpha = 1 - beta - if alpha[i, j, element] > new_alpha + newton_abstol - error("Alpha is getting smaller. old: $(alpha[i, j, element]), new: $new_alpha") - else - alpha[i, j, element] = new_alpha - end + alpha[i, j, element] = new_alpha return nothing end ### Auxiliary routines for Newton's bisection method ### # Initial checks -@inline function initial_check_entropy_spec_newton_idp(bound, goal, newton_abstol) +@inline function initial_check_local_onesided_newton_idp(::typeof(min), bound, + goal, newton_abstol) goal <= max(newton_abstol, abs(bound) * newton_abstol) end -@inline function initial_check_entropy_math_newton_idp(bound, goal, newton_abstol) +@inline function initial_check_local_onesided_newton_idp(::typeof(max), bound, + goal, newton_abstol) goal >= -max(newton_abstol, abs(bound) * newton_abstol) end -@inline initial_check_nonnegative_newton_idp(bound, goal, newton_abstol) = goal <= 0 +@inline initial_check_nonnegative_newton_idp(min_or_max, bound, goal, newton_abstol) = goal <= + 0 # Goal and d(Goal)d(u) function @inline goal_function_newton_idp(variable, bound, u, equations) = bound - @@ -671,8 +658,8 @@ end end # Final checks -# final check for entropy limiting -@inline function final_check_standard_newton_idp(bound, goal, newton_abstol) +# final check for one-sided local limiting +@inline function final_check_local_onesided_newton_idp(bound, goal, newton_abstol) abs(goal) < max(newton_abstol, abs(bound) * newton_abstol) end diff --git a/test/test_unit.jl b/test/test_unit.jl index ce735015012..716493a28c7 100644 --- a/test/test_unit.jl +++ b/test/test_unit.jl @@ -416,8 +416,9 @@ end indicator_hg = IndicatorHennemannGassner(1.0, 0.0, true, "variable", "cache") @test_nowarn show(stdout, indicator_hg) - limiter_idp = SubcellLimiterIDP(true, [1], true, [1], ["variable"], 0.1, true, true, - "cache", 1, (1.0, 1.0), 1.0) + limiter_idp = SubcellLimiterIDP(true, [1], true, [1], ["variable"], 0.1, + true, [(Trixi.entropy_spec, min)], "cache", + 1, (1.0, 1.0), 1.0) @test_nowarn show(stdout, limiter_idp) # TODO: TrixiShallowWater: move unit test From 6c7b2f596084380b150034a0aa8ad978578c3db6 Mon Sep 17 00:00:00 2001 From: bennibolm Date: Fri, 1 Mar 2024 14:36:16 +0100 Subject: [PATCH 16/25] Add comment to NEWS.md --- NEWS.md | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/NEWS.md b/NEWS.md index d70504d8c85..0c0fd8a8f7e 100644 --- a/NEWS.md +++ b/NEWS.md @@ -4,6 +4,18 @@ Trixi.jl follows the interpretation of [semantic versioning (semver)](https://ju used in the Julia ecosystem. Notable changes will be documented in this file for human readability. +## Changes in the v0.7 lifecycle + +#### Added +- Subcell local one-sided limiting support for nonlinear variables in 2D for `TreeMesh` + +#### Changed + +#### Deprecated + +#### Removed + + ## Changes when updating to v0.7 from v0.6.x #### Added From 6e5d9610bd029297a7673d7ad4851e440ca8fd60 Mon Sep 17 00:00:00 2001 From: bennibolm Date: Tue, 19 Mar 2024 14:34:31 +0100 Subject: [PATCH 17/25] Remove whitespaces --- NEWS.md | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/NEWS.md b/NEWS.md index 0ef5a594a96..00918fe0b43 100644 --- a/NEWS.md +++ b/NEWS.md @@ -17,7 +17,7 @@ for human readability. #### Changed -- The default wave speed estimate used within `flux_hll` is now `min_max_speed_davis` +- The default wave speed estimate used within `flux_hll` is now `min_max_speed_davis` instead of `min_max_speed_naive`. #### Deprecated @@ -33,7 +33,7 @@ for human readability. #### Added - AMR for hyperbolic-parabolic equations on 3D `P4estMesh` - `flux_hllc` on non-cartesian meshes for `CompressibleEulerEquations{2,3}D` -- Different boundary conditions for quad/hex meshes in Abaqus format, even if not generated by HOHQMesh, +- Different boundary conditions for quad/hex meshes in Abaqus format, even if not generated by HOHQMesh, can now be digested by Trixi in 2D and 3D. - Subcell (positivity) limiting support for nonlinear variables in 2D for `TreeMesh` - Added Lighthill-Whitham-Richards (LWR) traffic model @@ -47,7 +47,7 @@ for human readability. #### Changed - The wave speed estimates for `flux_hll`, `FluxHLL()` are now consistent across equations. - In particular, the functions `min_max_speed_naive`, `min_max_speed_einfeldt` are now + In particular, the functions `min_max_speed_naive`, `min_max_speed_einfeldt` are now conceptually identical across equations. Users, who have been using `flux_hll` for MHD have now to use `flux_hlle` in order to use the Einfeldt wave speed estimate. From ad7e1dabf3c532d18b8af45b0b490b47b4a4e3ef Mon Sep 17 00:00:00 2001 From: bennibolm Date: Mon, 25 Mar 2024 10:11:31 +0100 Subject: [PATCH 18/25] Implement suggestions --- src/callbacks_stage/subcell_bounds_check_2d.jl | 10 +++------- src/equations/compressible_euler_2d.jl | 15 +++------------ 2 files changed, 6 insertions(+), 19 deletions(-) diff --git a/src/callbacks_stage/subcell_bounds_check_2d.jl b/src/callbacks_stage/subcell_bounds_check_2d.jl index 4b8fb34c6b5..b8241bb8f72 100644 --- a/src/callbacks_stage/subcell_bounds_check_2d.jl +++ b/src/callbacks_stage/subcell_bounds_check_2d.jl @@ -47,18 +47,14 @@ foreach(limiter.local_onesided_variables_nonlinear) do (variable, min_or_max) key = Symbol(string(variable), "_", string(min_or_max)) deviation_threaded = idp_bounds_delta_local[key] + sign_ = min_or_max(1.0, -1.0) @threaded for element in eachelement(solver, cache) deviation = deviation_threaded[stride_size * Threads.threadid()] for j in eachnode(solver), i in eachnode(solver) v = variable(get_node_vars(u, equations, solver, i, j, element), equations) - if min_or_max === max - deviation = max(deviation, - v - variable_bounds[key][i, j, element]) - else # min_or_max === min - deviation = max(deviation, - variable_bounds[key][i, j, element] - v) - end + deviation = max(deviation, + sign_ * (v - variable_bounds[key][i, j, element])) end deviation_threaded[stride_size * Threads.threadid()] = deviation end diff --git a/src/equations/compressible_euler_2d.jl b/src/equations/compressible_euler_2d.jl index b6b83a2082f..9109bb24fa9 100644 --- a/src/equations/compressible_euler_2d.jl +++ b/src/equations/compressible_euler_2d.jl @@ -1854,8 +1854,9 @@ end return SVector(w1, w2, w3, w4) end -# Transformation from conservative variables u to entropy vector dSdu, -# S = -rho*s/(gamma-1), s=ln(p)-gamma*ln(rho) +# Transformation from conservative variables u to entropy vector ds_0/du, +# using the specific entropy of Guermond et al. (2019): s_0 = p * rho^(-gamma) / (gamma-1). +# Note: This is *not* the "conventional" specific entropy s = ln(p / rho^(gamma)). @inline function cons2entropy_spec(u, equations::CompressibleEulerEquations2D) rho, rho_v1, rho_v2, rho_e = u @@ -1871,13 +1872,6 @@ end w3 = -rho_v2 * inv_rho_gammap1 w4 = (1 / rho)^equations.gamma - # The derivative vector for other specific entropy - # sp = 1.0/(gammam1 * (rho_e - 0.5 * rho * v_square) - # w1 = gammam1 * 0.5 * v_square * sp - gamma / rho - # w2 = -gammam1 * v1 * sp - # w3 = -gammam1 * v2 * sp - # w4 = gammam1 * sp - return SVector(w1, w2, w3, w4) end @@ -1999,9 +1993,6 @@ end # Modified specific entropy from Guermond et al. (2019) s = (rho_e - 0.5 * (rho_v1^2 + rho_v2^2) / rho) * (1 / rho)^equations.gamma - # Other specific entropy - # rho_sp = rho/((equations.gamma - 1.0) * (rho_e - 0.5 * rho * v_square)) - # s = log(p) - (equaions.gamma + 1) * log(rho) return s end From 27a93f402d0cdeda25dd77d348e08b3add87c9a0 Mon Sep 17 00:00:00 2001 From: bennibolm Date: Mon, 25 Mar 2024 11:33:49 +0100 Subject: [PATCH 19/25] Replace `foreach` with `for` --- src/callbacks_stage/subcell_bounds_check_2d.jl | 2 +- src/solvers/dgsem_tree/subcell_limiters_2d.jl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/callbacks_stage/subcell_bounds_check_2d.jl b/src/callbacks_stage/subcell_bounds_check_2d.jl index b8241bb8f72..8e8cf025745 100644 --- a/src/callbacks_stage/subcell_bounds_check_2d.jl +++ b/src/callbacks_stage/subcell_bounds_check_2d.jl @@ -44,7 +44,7 @@ end end if local_onesided - foreach(limiter.local_onesided_variables_nonlinear) do (variable, min_or_max) + for (variable, min_or_max) in limiter.local_onesided_variables_nonlinear key = Symbol(string(variable), "_", string(min_or_max)) deviation_threaded = idp_bounds_delta_local[key] sign_ = min_or_max(1.0, -1.0) diff --git a/src/solvers/dgsem_tree/subcell_limiters_2d.jl b/src/solvers/dgsem_tree/subcell_limiters_2d.jl index e820bed667c..ac75b88de29 100644 --- a/src/solvers/dgsem_tree/subcell_limiters_2d.jl +++ b/src/solvers/dgsem_tree/subcell_limiters_2d.jl @@ -354,7 +354,7 @@ end # Local one-sided limiting of nonlinear variables @inline function idp_local_onesided!(alpha, limiter, u, t, dt, semi) - foreach(limiter.local_onesided_variables_nonlinear) do (variable, min_or_max) + for (variable, min_or_max) in limiter.local_onesided_variables_nonlinear idp_local_onesided!(alpha, limiter, u, t, dt, semi, variable, min_or_max) end From 673211b7787e07e2fab079629b58f8485f38612e Mon Sep 17 00:00:00 2001 From: bennibolm Date: Mon, 25 Mar 2024 14:38:46 +0100 Subject: [PATCH 20/25] Fix tests --- test/test_tree_2d_euler.jl | 32 ++++++++++++++++---------------- 1 file changed, 16 insertions(+), 16 deletions(-) diff --git a/test/test_tree_2d_euler.jl b/test/test_tree_2d_euler.jl index 44feab9f40f..87898953a71 100644 --- a/test/test_tree_2d_euler.jl +++ b/test/test_tree_2d_euler.jl @@ -317,16 +317,16 @@ end @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_blast_wave_sc_subcell_nonperiodic.jl"), l2=[ - 0.32211935909459755, - 0.17984734129031393, - 0.17983517637561672, - 0.6136858672602447, + 0.322114786967044, + 0.17984537743598789, + 0.17983274405937094, + 0.6136851154934142, ], linf=[ - 1.3435940880493509, - 1.1748248970276045, - 1.1745481442875036, - 2.421529617190895, + 1.3435381344401196, + 1.1748892702276696, + 1.1744575084504332, + 2.4215286889351857, ], tspan=(0.0, 0.5), initial_refinement_level=4, @@ -371,16 +371,16 @@ end @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_sedov_blast_wave_sc_subcell.jl"), l2=[ - 0.4144617128623457, - 0.14606474605745345, - 0.14607706120162497, - 0.616836505667339, + 0.41446170615771283, + 0.1460647720630636, + 0.14607708378676684, + 0.61683646759431, ], linf=[ - 1.5717812331790788, - 0.7963178376801539, - 0.795685946979518, - 6.456116356402219, + 1.5717833640175574, + 0.796318274890205, + 0.7956868775911277, + 6.456117275380676, ], tspan=(0.0, 1.0), initial_refinement_level=4, From f059814770f89da5ecdfde2710f5e6b6d6a4252d Mon Sep 17 00:00:00 2001 From: bennibolm Date: Mon, 8 Apr 2024 11:12:46 +0200 Subject: [PATCH 21/25] Use new bounds check reduction for one sided limiter --- src/callbacks_stage/subcell_bounds_check_2d.jl | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/src/callbacks_stage/subcell_bounds_check_2d.jl b/src/callbacks_stage/subcell_bounds_check_2d.jl index f1ddadd2a77..0f713a296e2 100644 --- a/src/callbacks_stage/subcell_bounds_check_2d.jl +++ b/src/callbacks_stage/subcell_bounds_check_2d.jl @@ -48,18 +48,20 @@ if local_onesided for (variable, min_or_max) in limiter.local_onesided_variables_nonlinear key = Symbol(string(variable), "_", string(min_or_max)) - deviation_threaded = idp_bounds_delta_local[key] + deviation = idp_bounds_delta_local[key] sign_ = min_or_max(1.0, -1.0) - @threaded for element in eachelement(solver, cache) - deviation = deviation_threaded[stride_size * Threads.threadid()] + @batch reduction=(max, deviation) for element in eachelement(solver, cache) for j in eachnode(solver), i in eachnode(solver) v = variable(get_node_vars(u, equations, solver, i, j, element), equations) + # Note: We always save the absolute deviations >= 0 and therefore use the + # `max` operator for lower and upper bounds. The different directions of + # upper and lower bounds are considered with `sign_`. deviation = max(deviation, sign_ * (v - variable_bounds[key][i, j, element])) end - deviation_threaded[stride_size * Threads.threadid()] = deviation end + idp_bounds_delta_local[key] = deviation end end if positivity From da1af1e29d2134c8a96b0595e7b5875a40c35b7d Mon Sep 17 00:00:00 2001 From: bennibolm Date: Tue, 7 May 2024 12:11:15 +0200 Subject: [PATCH 22/25] Adapt tests after last merge of main --- test/test_tree_2d_euler.jl | 32 ++++++++++++++++---------------- 1 file changed, 16 insertions(+), 16 deletions(-) diff --git a/test/test_tree_2d_euler.jl b/test/test_tree_2d_euler.jl index 2c82eea4a5a..efd4058031f 100644 --- a/test/test_tree_2d_euler.jl +++ b/test/test_tree_2d_euler.jl @@ -349,16 +349,16 @@ end @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_blast_wave_sc_subcell_nonperiodic.jl"), l2=[ - 0.322114786967044, - 0.17984537743598789, - 0.17983274405937094, - 0.6136851154934142, + 0.3221177942225801, + 0.1798478357478982, + 0.1798364616438908, + 0.6136884131056267, ], linf=[ - 1.3435381344401196, - 1.1748892702276696, - 1.1744575084504332, - 2.4215286889351857, + 1.343766644801395, + 1.1749593109683463, + 1.1747613085307178, + 2.4216006041018785, ], tspan=(0.0, 0.5), initial_refinement_level=4, @@ -403,16 +403,16 @@ end @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_sedov_blast_wave_sc_subcell.jl"), l2=[ - 0.41446170615771283, - 0.1460647720630636, - 0.14607708378676684, - 0.61683646759431, + 0.41444427153173785, + 0.1460669409661223, + 0.14606693069201596, + 0.6168046457461059, ], linf=[ - 1.5717833640175574, - 0.796318274890205, - 0.7956868775911277, - 6.456117275380676, + 1.5720584643579567, + 0.7946656826861964, + 0.7946656525739751, + 6.455520291414711, ], tspan=(0.0, 1.0), initial_refinement_level=4, From 27ded75f469a66d6632815626d373549710bf1fc Mon Sep 17 00:00:00 2001 From: bennibolm Date: Tue, 7 May 2024 14:19:18 +0200 Subject: [PATCH 23/25] Rename `entropy_spec` to `entropy_guermond_etal` --- .../elixir_euler_sedov_blast_wave_sc_subcell.jl | 2 +- src/equations/compressible_euler_2d.jl | 13 +++++++------ src/solvers/dgsem_tree/subcell_limiters.jl | 6 +++--- test/test_unit.jl | 2 +- 4 files changed, 12 insertions(+), 11 deletions(-) diff --git a/examples/tree_2d_dgsem/elixir_euler_sedov_blast_wave_sc_subcell.jl b/examples/tree_2d_dgsem/elixir_euler_sedov_blast_wave_sc_subcell.jl index d3bb6a819b5..6cbbe4eb4e6 100644 --- a/examples/tree_2d_dgsem/elixir_euler_sedov_blast_wave_sc_subcell.jl +++ b/examples/tree_2d_dgsem/elixir_euler_sedov_blast_wave_sc_subcell.jl @@ -43,7 +43,7 @@ volume_flux = flux_chandrashekar basis = LobattoLegendreBasis(3) limiter_idp = SubcellLimiterIDP(equations, basis; local_twosided_variables_cons = ["rho"], - local_onesided_variables_nonlinear = [(Trixi.entropy_spec, + local_onesided_variables_nonlinear = [(Trixi.entropy_guermond_etal, min)]) volume_integral = VolumeIntegralSubcellLimiting(limiter_idp; volume_flux_dg = volume_flux, diff --git a/src/equations/compressible_euler_2d.jl b/src/equations/compressible_euler_2d.jl index 3cc94415d70..0614066806c 100644 --- a/src/equations/compressible_euler_2d.jl +++ b/src/equations/compressible_euler_2d.jl @@ -1887,9 +1887,9 @@ end end # Transformation from conservative variables u to entropy vector ds_0/du, -# using the specific entropy of Guermond et al. (2019): s_0 = p * rho^(-gamma) / (gamma-1). +# using the modified specific entropy of Guermond et al. (2019): s_0 = p * rho^(-gamma) / (gamma-1). # Note: This is *not* the "conventional" specific entropy s = ln(p / rho^(gamma)). -@inline function cons2entropy_spec(u, equations::CompressibleEulerEquations2D) +@inline function cons2entropy_guermond_etal(u, equations::CompressibleEulerEquations2D) rho, rho_v1, rho_v2, rho_e = u v1 = rho_v1 / rho @@ -2018,8 +2018,9 @@ end return cons2entropy(u, equations) end -# Calculate specific entropy for conservative variable u -@inline function entropy_spec(u, equations::CompressibleEulerEquations2D) +# Calculate the modified specific entropy of Guermond et al. (2019): s_0 = p * rho^(-gamma) / (gamma-1). +# Note: This is *not* the "conventional" specific entropy s = ln(p / rho^(gamma)). +@inline function entropy_guermond_etal(u, equations::CompressibleEulerEquations2D) rho, rho_v1, rho_v2, rho_e = u # Modified specific entropy from Guermond et al. (2019) @@ -2029,9 +2030,9 @@ end end # Transformation from conservative variables u to d(s)/d(u) -@inline function gradient_conservative(::typeof(entropy_spec), +@inline function gradient_conservative(::typeof(entropy_guermond_etal), u, equations::CompressibleEulerEquations2D) - return cons2entropy_spec(u, equations) + return cons2entropy_guermond_etal(u, equations) end # Default entropy is the mathematical entropy diff --git a/src/solvers/dgsem_tree/subcell_limiters.jl b/src/solvers/dgsem_tree/subcell_limiters.jl index 7b06149148d..d29c6f2e39b 100644 --- a/src/solvers/dgsem_tree/subcell_limiters.jl +++ b/src/solvers/dgsem_tree/subcell_limiters.jl @@ -28,7 +28,7 @@ including: - Local two-sided Zalesak-type limiting for conservative variables (`local_twosided_variables_cons`) - Positivity limiting for conservative variables (`positivity_variables_cons`) and nonlinear variables (`positivity_variables_nonlinear`) -- Local one-sided limiting for nonlinear variables, e.g. `entropy_spec` and `entropy_math` +- Local one-sided limiting for nonlinear variables, e.g. `entropy_guermond_etal` and `entropy_math` with `local_onesided_variables_nonlinear` To use these three limiting options use the following structure: @@ -38,8 +38,8 @@ To use these three limiting options use the following structure: For ***nonlinear variables***, the wanted variable functions are passed within a vector: To ensure positivity use a plain vector including the desired variables, e.g. `positivity_variables_nonlinear = [pressure]`. For local one-sided limiting pass the variable function combined with the requested bound -(`min` or `max`) as a tuple. For instance, to impose a lower local bound on the specific entropy -use `local_onesided_variables_nonlinear = [(Trixi.entropy_spec, min)]`. +(`min` or `max`) as a tuple. For instance, to impose a lower local bound on the modified specific +entropy by Guermond et al. use `local_onesided_variables_nonlinear = [(Trixi.entropy_guermond_etal, min)]`. The bounds are calculated using the low-order FV solution. The positivity limiter uses `positivity_correction_factor` such that `u^new >= positivity_correction_factor * u^FV`. diff --git a/test/test_unit.jl b/test/test_unit.jl index f9ef07384c4..90ee21030d3 100644 --- a/test/test_unit.jl +++ b/test/test_unit.jl @@ -417,7 +417,7 @@ end @test_nowarn show(stdout, indicator_hg) limiter_idp = SubcellLimiterIDP(true, [1], true, [1], ["variable"], 0.1, - true, [(Trixi.entropy_spec, min)], "cache", + true, [(Trixi.entropy_guermond_etal, min)], "cache", 1, (1.0, 1.0), 1.0) @test_nowarn show(stdout, limiter_idp) From 62946c2220b1005ac1ff9e8976425590d22e27e9 Mon Sep 17 00:00:00 2001 From: bennibolm Date: Tue, 7 May 2024 14:19:47 +0200 Subject: [PATCH 24/25] Remove not-needed `tuple` --- src/solvers/dgsem_tree/subcell_limiters.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/solvers/dgsem_tree/subcell_limiters.jl b/src/solvers/dgsem_tree/subcell_limiters.jl index d29c6f2e39b..17c9488316d 100644 --- a/src/solvers/dgsem_tree/subcell_limiters.jl +++ b/src/solvers/dgsem_tree/subcell_limiters.jl @@ -101,11 +101,11 @@ function SubcellLimiterIDP(equations::AbstractEquations, basis; local_onesided_variables_nonlinear_ = Tuple{Function, Function}[] for (variable, min_or_max) in local_onesided_variables_nonlinear if min_or_max === Base.max - push!(local_onesided_variables_nonlinear_, tuple(variable, max)) + push!(local_onesided_variables_nonlinear_, (variable, max)) elseif min_or_max === Base.min - push!(local_onesided_variables_nonlinear_, tuple(variable, min)) + push!(local_onesided_variables_nonlinear_, (variable, min)) elseif min_or_max === Trixi.max || min_or_max === Trixi.min - push!(local_onesided_variables_nonlinear_, tuple(variable, min_or_max)) + push!(local_onesided_variables_nonlinear_, (variable, min_or_max)) else error("Parameter $min_or_max is not a valid input. Use `max` or `min` instead.") end From 03f9f0545ec13532e4465556a800a194a4ac78a7 Mon Sep 17 00:00:00 2001 From: bennibolm Date: Tue, 7 May 2024 14:54:07 +0200 Subject: [PATCH 25/25] Adapt nameing changes to tutorial --- docs/literate/src/files/subcell_shock_capturing.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/docs/literate/src/files/subcell_shock_capturing.jl b/docs/literate/src/files/subcell_shock_capturing.jl index 8b8e49a28a8..8a98fdae283 100644 --- a/docs/literate/src/files/subcell_shock_capturing.jl +++ b/docs/literate/src/files/subcell_shock_capturing.jl @@ -106,7 +106,7 @@ positivity_variables_nonlinear = [pressure] # As for the limiting with global bounds you are passing the quantity names of the conservative # variables you want to limit. So, to limit the density with lower and upper local bounds pass # the following. -local_minmax_variables_cons = ["rho"] +local_twosided_variables_cons = ["rho"] # ## Exemplary simulation # How to set up a simulation using the IDP limiting becomes clearer when looking at an exemplary @@ -154,7 +154,7 @@ volume_flux = flux_ranocha # Here, the simulation should contain local limiting for the density using lower and upper bounds. basis = LobattoLegendreBasis(3) limiter_idp = SubcellLimiterIDP(equations, basis; - local_minmax_variables_cons = ["rho"]) + local_twosided_variables_cons = ["rho"]) # The initialized limiter is passed to `VolumeIntegralSubcellLimiting` in addition to the volume # fluxes of the low-order and high-order scheme.