From 7acf81b80b139a83212347e72e000e5cdafabf46 Mon Sep 17 00:00:00 2001 From: bennibolm Date: Thu, 13 Nov 2025 11:33:35 +0100 Subject: [PATCH 1/5] Add 3d positivity limiting for nonlinear variables --- .../elixir_euler_sedov_sc_subcell.jl | 7 +- src/equations/compressible_euler_3d.jl | 9 ++ .../dgsem_p4est/subcell_limiters_3d.jl | 103 +++++++++++++++++- test/test_p4est_3d.jl | 2 +- 4 files changed, 116 insertions(+), 5 deletions(-) diff --git a/examples/p4est_3d_dgsem/elixir_euler_sedov_sc_subcell.jl b/examples/p4est_3d_dgsem/elixir_euler_sedov_sc_subcell.jl index b8be1d7d85c..e9aa93ecc34 100644 --- a/examples/p4est_3d_dgsem/elixir_euler_sedov_sc_subcell.jl +++ b/examples/p4est_3d_dgsem/elixir_euler_sedov_sc_subcell.jl @@ -25,7 +25,8 @@ function initial_condition_sedov_blast_wave(x, t, r0 = 0.21875 # = 3.5 * smallest dx (for domain length=4 and max-ref=6) E = 1.0 p0_inner = 3 * (equations.gamma - 1) * E / (4 * pi * r0^2) - p0_outer = 1.0e-1 # "simpler" setup since positivity limiter for pressure is not yet supported in 3D + p0_outer = 1.0e-5 # = true Sedov setup + # p0_outer = 1.0e-3 # = more reasonable setup # Calculate primitive variables rho = 1.0 @@ -36,7 +37,6 @@ function initial_condition_sedov_blast_wave(x, t, return prim2cons(SVector(rho, v1, v2, v3, p), equations) end - initial_condition = initial_condition_sedov_blast_wave surface_flux = flux_lax_friedrichs @@ -44,7 +44,8 @@ volume_flux = flux_ranocha polydeg = 3 basis = LobattoLegendreBasis(polydeg) limiter_idp = SubcellLimiterIDP(equations, basis; - positivity_variables_cons = ["rho"]) + positivity_variables_cons = ["rho"], + positivity_variables_nonlinear = [pressure]) volume_integral = VolumeIntegralSubcellLimiting(limiter_idp; volume_flux_dg = volume_flux, volume_flux_fv = surface_flux) diff --git a/src/equations/compressible_euler_3d.jl b/src/equations/compressible_euler_3d.jl index 41a93001c9d..afedf924974 100644 --- a/src/equations/compressible_euler_3d.jl +++ b/src/equations/compressible_euler_3d.jl @@ -1898,4 +1898,13 @@ end @inline function energy_internal(cons, equations::CompressibleEulerEquations3D) return energy_total(cons, equations) - energy_kinetic(cons, equations) end + +# State validation for Newton-bisection method of subcell IDP limiting +@inline function Base.isvalid(u, equations::CompressibleEulerEquations3D) + p = pressure(u, equations) + if u[1] <= 0 || p <= 0 + return false + end + return true +end end # @muladd diff --git a/src/solvers/dgsem_p4est/subcell_limiters_3d.jl b/src/solvers/dgsem_p4est/subcell_limiters_3d.jl index 492e85ab188..8a60d349b88 100644 --- a/src/solvers/dgsem_p4est/subcell_limiters_3d.jl +++ b/src/solvers/dgsem_p4est/subcell_limiters_3d.jl @@ -29,7 +29,7 @@ mesh, i, j, k, element) var = u[variable, i, j, k, element] if var < 0 - error("Safe low-order method produces negative value for conservative variable $variable. Try a smaller time step.") + throw(ArgumentError("Safe low-order method produces negative value for conservative variable $variable. Try a smaller time step.")) end # Compute bound @@ -73,4 +73,105 @@ return nothing end + +############################################################################### +# Global positivity limiting of nonlinear variables + +@inline function idp_positivity_nonlinear!(alpha, limiter, + u::AbstractArray{<:Real, 5}, 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 k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) + inverse_jacobian = get_inverse_jacobian(cache.elements.inverse_jacobian, + mesh, i, j, k, element) + + # Compute bound + u_local = get_node_vars(u, equations, dg, i, j, k, element) + var = variable(u_local, equations) + if var < 0 + throw(ArgumentError("Safe low-order method produces negative value for variable $variable. Try a smaller time step.")) + end + var_min[i, j, k, element] = positivity_correction_factor * var + + # Perform Newton's bisection method to find new alpha + newton_loops_alpha!(alpha, var_min[i, j, k, element], + u_local, i, j, k, element, + variable, min, + initial_check_nonnegative_newton_idp, + final_check_nonnegative_newton_idp, + inverse_jacobian, dt, equations, dg, cache, limiter) + end + end + + return nothing +end + +############################################################################### +# Newton-bisection method + +@inline function newton_loops_alpha!(alpha, bound, u, i, j, k, element, + variable, min_or_max, + initial_check, final_check, + inverse_jacobian, dt, + equations, dg, cache, limiter) + (; inverse_weights) = dg.basis + (; antidiffusive_flux1_L, antidiffusive_flux1_R, antidiffusive_flux2_L, antidiffusive_flux2_R, antidiffusive_flux3_L, antidiffusive_flux3_R) = cache.antidiffusive_fluxes + + (; gamma_constant_newton) = limiter + + indices = (i, j, k, element) + + # negative xi direction + antidiffusive_flux = gamma_constant_newton * inverse_jacobian * inverse_weights[i] * + get_node_vars(antidiffusive_flux1_R, equations, dg, + i, j, k, element) + newton_loop!(alpha, bound, u, indices, 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, k, element) + newton_loop!(alpha, bound, u, indices, 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, k, element) + newton_loop!(alpha, bound, u, indices, 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, k, element) + newton_loop!(alpha, bound, u, indices, variable, min_or_max, + initial_check, final_check, equations, dt, limiter, antidiffusive_flux) + + # negative zeta direction + antidiffusive_flux = gamma_constant_newton * inverse_jacobian * inverse_weights[k] * + get_node_vars(antidiffusive_flux3_R, equations, dg, + i, j, k, element) + newton_loop!(alpha, bound, u, indices, variable, min_or_max, + initial_check, final_check, equations, dt, limiter, antidiffusive_flux) + + # positive zeta direction + antidiffusive_flux = -gamma_constant_newton * inverse_jacobian * + inverse_weights[k] * + get_node_vars(antidiffusive_flux3_L, equations, dg, + i, j, k + 1, element) + newton_loop!(alpha, bound, u, indices, variable, min_or_max, + initial_check, final_check, equations, dt, limiter, antidiffusive_flux) + + return nothing +end end # @muladd diff --git a/test/test_p4est_3d.jl b/test/test_p4est_3d.jl index 42aba382f89..e3f5fba22b7 100644 --- a/test/test_p4est_3d.jl +++ b/test/test_p4est_3d.jl @@ -340,7 +340,7 @@ end 2.8814684266647856, 6.002617193075404 ], - tspan=(0.0, 0.2),) + tspan=(0.0, 0.3),) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) # Larger values for allowed allocations due to usage of custom From 2bb91dd9d3d93cc614a01381556d570c305aff14 Mon Sep 17 00:00:00 2001 From: bennibolm Date: Thu, 13 Nov 2025 12:05:52 +0100 Subject: [PATCH 2/5] Adapt test results --- test/test_p4est_3d.jl | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/test/test_p4est_3d.jl b/test/test_p4est_3d.jl index e3f5fba22b7..21ee330435b 100644 --- a/test/test_p4est_3d.jl +++ b/test/test_p4est_3d.jl @@ -327,18 +327,18 @@ end @trixi_testset "elixir_euler_sedov_sc_subcell.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_sedov_sc_subcell.jl"), l2=[ - 0.16371327294681623, - 0.0779033135575244, - 0.07790331355752474, - 0.07790331355752393, - 0.36689870760596666 + 0.1942700455652903, + 0.07557644365785855, + 0.07557644365785836, + 0.07557644365785698, + 0.3713893635249306 ], linf=[ - 4.210879608740924, - 2.881468426664787, - 2.8814684266647874, - 2.8814684266647856, - 6.002617193075404 + 2.7542157588958798, + 1.8885700263691245, + 1.888570026369125, + 1.8885700263691252, + 4.9712792944452096 ], tspan=(0.0, 0.3),) # Ensure that we do not have excessive memory allocations From a2bccd38c24b86b3e3e1b919bbc5bae9a3662b53 Mon Sep 17 00:00:00 2001 From: bennibolm Date: Mon, 17 Nov 2025 09:15:37 +0100 Subject: [PATCH 3/5] Implement suggestions --- examples/p4est_3d_dgsem/elixir_euler_sedov_sc_subcell.jl | 1 - src/equations/compressible_euler_2d.jl | 3 +-- src/equations/compressible_euler_3d.jl | 3 +-- src/solvers/dgsem_p4est/subcell_limiters_3d.jl | 6 +++--- 4 files changed, 5 insertions(+), 8 deletions(-) diff --git a/examples/p4est_3d_dgsem/elixir_euler_sedov_sc_subcell.jl b/examples/p4est_3d_dgsem/elixir_euler_sedov_sc_subcell.jl index e9aa93ecc34..62186cd1792 100644 --- a/examples/p4est_3d_dgsem/elixir_euler_sedov_sc_subcell.jl +++ b/examples/p4est_3d_dgsem/elixir_euler_sedov_sc_subcell.jl @@ -26,7 +26,6 @@ function initial_condition_sedov_blast_wave(x, t, E = 1.0 p0_inner = 3 * (equations.gamma - 1) * E / (4 * pi * r0^2) p0_outer = 1.0e-5 # = true Sedov setup - # p0_outer = 1.0e-3 # = more reasonable setup # Calculate primitive variables rho = 1.0 diff --git a/src/equations/compressible_euler_2d.jl b/src/equations/compressible_euler_2d.jl index 81b63781e54..2c8341f9a61 100644 --- a/src/equations/compressible_euler_2d.jl +++ b/src/equations/compressible_euler_2d.jl @@ -2232,8 +2232,7 @@ end # State validation for Newton-bisection method of subcell IDP limiting @inline function Base.isvalid(u, equations::CompressibleEulerEquations2D) - p = pressure(u, equations) - if u[1] <= 0 || p <= 0 + if u[1] <= 0 || pressure(u, equations) <= 0 return false end return true diff --git a/src/equations/compressible_euler_3d.jl b/src/equations/compressible_euler_3d.jl index afedf924974..4258735c8d5 100644 --- a/src/equations/compressible_euler_3d.jl +++ b/src/equations/compressible_euler_3d.jl @@ -1901,8 +1901,7 @@ end # State validation for Newton-bisection method of subcell IDP limiting @inline function Base.isvalid(u, equations::CompressibleEulerEquations3D) - p = pressure(u, equations) - if u[1] <= 0 || p <= 0 + if u[1] <= 0 || pressure(u, equations) <= 0 return false end return true diff --git a/src/solvers/dgsem_p4est/subcell_limiters_3d.jl b/src/solvers/dgsem_p4est/subcell_limiters_3d.jl index 8a60d349b88..4a83e68f6f0 100644 --- a/src/solvers/dgsem_p4est/subcell_limiters_3d.jl +++ b/src/solvers/dgsem_p4est/subcell_limiters_3d.jl @@ -29,7 +29,7 @@ mesh, i, j, k, element) var = u[variable, i, j, k, element] if var < 0 - throw(ArgumentError("Safe low-order method produces negative value for conservative variable $variable. Try a smaller time step.")) + error("Safe low-order method produces negative value for conservative variable $variable. Try a smaller time step.") end # Compute bound @@ -95,7 +95,7 @@ end u_local = get_node_vars(u, equations, dg, i, j, k, element) var = variable(u_local, equations) if var < 0 - throw(ArgumentError("Safe low-order method produces negative value for variable $variable. Try a smaller time step.")) + error("Safe low-order method produces negative value for variable $variable. Try a smaller time step.") end var_min[i, j, k, element] = positivity_correction_factor * var @@ -120,7 +120,7 @@ end initial_check, final_check, inverse_jacobian, dt, equations, dg, cache, limiter) - (; inverse_weights) = dg.basis + (; inverse_weights) = dg.basis # Plays role of inverse DG-subcell sizes (; antidiffusive_flux1_L, antidiffusive_flux1_R, antidiffusive_flux2_L, antidiffusive_flux2_R, antidiffusive_flux3_L, antidiffusive_flux3_R) = cache.antidiffusive_fluxes (; gamma_constant_newton) = limiter From 0e96e129d29ad959c2f2d0eaecda23f640d901f8 Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Mon, 17 Nov 2025 11:26:16 +0100 Subject: [PATCH 4/5] Apply suggestions from code review --- src/solvers/dgsem_p4est/subcell_limiters_3d.jl | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/src/solvers/dgsem_p4est/subcell_limiters_3d.jl b/src/solvers/dgsem_p4est/subcell_limiters_3d.jl index 4a83e68f6f0..b5f6b729f26 100644 --- a/src/solvers/dgsem_p4est/subcell_limiters_3d.jl +++ b/src/solvers/dgsem_p4est/subcell_limiters_3d.jl @@ -128,7 +128,8 @@ end indices = (i, j, k, element) # negative xi direction - antidiffusive_flux = gamma_constant_newton * inverse_jacobian * inverse_weights[i] * + antidiffusive_flux = gamma_constant_newton * inverse_jacobian * + inverse_weights[i] * get_node_vars(antidiffusive_flux1_R, equations, dg, i, j, k, element) newton_loop!(alpha, bound, u, indices, variable, min_or_max, @@ -143,7 +144,8 @@ end initial_check, final_check, equations, dt, limiter, antidiffusive_flux) # negative eta direction - antidiffusive_flux = gamma_constant_newton * inverse_jacobian * inverse_weights[j] * + antidiffusive_flux = gamma_constant_newton * inverse_jacobian * + inverse_weights[j] * get_node_vars(antidiffusive_flux2_R, equations, dg, i, j, k, element) newton_loop!(alpha, bound, u, indices, variable, min_or_max, @@ -158,7 +160,8 @@ end initial_check, final_check, equations, dt, limiter, antidiffusive_flux) # negative zeta direction - antidiffusive_flux = gamma_constant_newton * inverse_jacobian * inverse_weights[k] * + antidiffusive_flux = gamma_constant_newton * inverse_jacobian * + inverse_weights[k] * get_node_vars(antidiffusive_flux3_R, equations, dg, i, j, k, element) newton_loop!(alpha, bound, u, indices, variable, min_or_max, From f7b4e003ff6a721bd4e727edf0c5f585fdd817d9 Mon Sep 17 00:00:00 2001 From: bennibolm Date: Wed, 19 Nov 2025 17:25:38 +0100 Subject: [PATCH 5/5] Update NEWS.md --- NEWS.md | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/NEWS.md b/NEWS.md index 5fe40e12792..c4a886cbf88 100644 --- a/NEWS.md +++ b/NEWS.md @@ -9,9 +9,10 @@ for human readability. ## Changes in the v0.13 lifecycle #### Added -- Initial 3D support for subcell limiting with `P4estMesh` was added ([#2582]). +- Initial 3D support for subcell limiting with `P4estMesh` was added ([#2582] and [#2647]). In the new version, IDP positivity limiting for conservative variables (using - the keyword `positivity_variables_cons` in `SubcellLimiterIDP()`) is supported. + the keyword `positivity_variables_cons` in `SubcellLimiterIDP()`) and nonlinear + variables (using `positivity_variables_nonlinear`) is supported. `BoundsCheckCallback` is not supported in 3D yet.