Add 3d subcell positivity limiting for nonlinear variables

examples/p4est_3d_dgsem/elixir_euler_sedov_sc_subcell.jl

@@ -25,7 +25,7 @@ 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
 
     # Calculate primitive variables
     rho = 1.0
@@ -36,15 +36,15 @@ 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
 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)

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

src/equations/compressible_euler_3d.jl

@@ -1898,4 +1898,12 @@ 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)
+    if u[1] <= 0 || pressure(u, equations) <= 0
+        return false
+    end
+    return true
+end
 end # @muladd

src/solvers/dgsem_p4est/subcell_limiters_3d.jl

@@ -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
+                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
+
+            # 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 # 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
+
+    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

test/test_p4est_3d.jl

@@ -327,20 +327,20 @@ 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.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