Add 3d local subcell limiting for nonlinear variables

NEWS.md

@@ -13,6 +13,9 @@ for human readability.
 - Support for second-order finite volume subcell volume integral (`VolumeIntegralPureLGLFiniteVolumeO2`) and
   stabilized DG-FV blending volume integral (`VolumeIntegralShockCapturingRRG`) on 
   3D meshes for conservative systems ([#2734], [#2755]).
+- Extended 3D support for subcell limiting with `P4estMesh` was added ([#2733]).
+  In the new version, local (minimum or maximum) limiting for nonlinear variables (using
+  the keyword `local_onesided_variables_nonlinear` in `SubcellLimiterIDP()`) is supported.
 
 #### Changed
 
@@ -26,17 +29,16 @@ for human readability.
 ## Changes in the v0.13 lifecycle
 
 #### Added
-- Initial 3D support for subcell limiting with `P4estMesh` was added ([#2582] and [#2647]).
+- Initial 3D support for subcell limiting with `P4estMesh` was added ([#2582], [#2647], [#2688], [#2722]).
   In the new version, IDP positivity limiting for conservative variables (using
   the keyword `positivity_variables_cons` in `SubcellLimiterIDP()`) and nonlinear
   variables (using `positivity_variables_nonlinear`) is supported.
-  `BoundsCheckCallback` is not supported in 3D yet.
 - Optimized 2D and 3D kernels for nonconservative fluxes with `P4estMesh` were added ([#2653], [#2663]).
   The optimized kernel can be enabled via the trait `Trixi.combine_conservative_and_nonconservative_fluxes(flux, equations)`.
   When the trait is set to `Trixi.True()`, a single method has to be defined, that computes and returns the tuple
   `flux_cons(u_ll, u_rr) + 0.5f0 * flux_noncons(u_ll, u_rr)` and
   `flux_cons(u_ll, u_rr) + 0.5f0 * flux_noncons(u_rr, u_ll)`.
-- Added support for `source_terms_parabolic`, which allows users to specify gradient-dependent source terms when solving parabolic equations ([#2721]). 
+- Added support for `source_terms_parabolic`, which allows users to specify gradient-dependent source terms when solving parabolic equations ([#2721]).
 - Support for second-order finite volume subcell volume integral (`VolumeIntegralPureLGLFiniteVolumeO2`) and
   stabilized DG-FV blending volume integral (`VolumeIntegralShockCapturingRRG`) on 
   1D & 2D meshes for conservative systems ([#2022], [#2695], [#2718]).

examples/p4est_3d_dgsem/elixir_euler_sedov_sc_subcell.jl

@@ -44,7 +44,9 @@ polydeg = 3
 basis = LobattoLegendreBasis(polydeg)
 limiter_idp = SubcellLimiterIDP(equations, basis;
                                 positivity_variables_cons = ["rho"],
-                                positivity_variables_nonlinear = [pressure])
+                                positivity_variables_nonlinear = [pressure],
+                                local_onesided_variables_nonlinear = [],
+                                max_iterations_newton = 25)
 volume_integral = VolumeIntegralSubcellLimiting(limiter_idp;
                                                 volume_flux_dg = volume_flux,
                                                 volume_flux_fv = surface_flux)

examples/p4est_3d_dgsem/elixir_mhd_shockcapturing_subcell.jl

@@ -117,7 +117,7 @@ callbacks = CallbackSet(summary_callback,
 
 ###############################################################################
 # run the simulation
-stage_callbacks = (SubcellLimiterIDPCorrection(),)
+stage_callbacks = (SubcellLimiterIDPCorrection(), BoundsCheckCallback())
 
 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

src/callbacks_stage/subcell_bounds_check_3d.jl

@@ -19,6 +19,26 @@
     # `@batch` here to allow a possible redefinition of `@threaded` without creating errors here.
     # See also https://github.com/trixi-framework/Trixi.jl/pull/1888#discussion_r1537785293.
 
+    if local_onesided
+        for (variable, min_or_max) in limiter.local_onesided_variables_nonlinear
+            key = Symbol(string(variable), "_", string(min_or_max))
+            deviation = idp_bounds_delta_local[key]
+            sign_ = min_or_max(1.0, -1.0)
+            @batch reduction=(max, deviation) for element in eachelement(solver, cache)
+                for k in eachnode(solver), j in eachnode(solver), i in eachnode(solver)
+                    v = variable(get_node_vars(u, equations, solver, i, j, k, 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, k, element]))
+                end
+            end
+            idp_bounds_delta_local[key] = deviation
+        end
+    end
     if positivity
         for v in limiter.positivity_variables_cons
             key = Symbol(string(v), "_min")

src/equations/compressible_euler_3d.jl

@@ -1791,6 +1791,27 @@ end
     return SVector(w1, w2, w3, w4, w5)
 end
 
+# Transformation from conservative variables u to entropy vector ds_0/du,
+# 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_guermond_etal(u, equations::CompressibleEulerEquations3D)
+    rho, rho_v1, rho_v2, rho_v3, rho_e = u
+
+    inv_rho = 1 / rho
+    v_square_rho = (rho_v1^2 + rho_v2^2 + rho_v3^2) * inv_rho
+    inv_rho_gammap1 = inv_rho^(equations.gamma + 1)
+
+    # The derivative vector for the modified specific entropy of Guermond et al.
+    w1 = inv_rho_gammap1 *
+         (0.5f0 * (equations.gamma + 1) * v_square_rho - equations.gamma * rho_e)
+    w2 = -rho_v1 * inv_rho_gammap1
+    w3 = -rho_v2 * inv_rho_gammap1
+    w4 = -rho_v3 * inv_rho_gammap1
+    w5 = inv_rho^equations.gamma
+
+    return SVector(w1, w2, w3, w4, w5)
+end
+
 @inline function entropy2cons(w, equations::CompressibleEulerEquations3D)
     # 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)
@@ -1880,6 +1901,34 @@ end
     return S
 end
 
+@doc raw"""
+    entropy_guermond_etal(u, equations::CompressibleEulerEquations3D)
+
+Calculate the modified specific entropy of Guermond et al. (2019):
+```math
+s_0 = p * \rho^{-\gamma} / (\gamma-1).
+```
+Note: This is *not* the "conventional" specific entropy ``s = ln(p / \rho^\gamma)``.
+- Guermond at al. (2019)
+  Invariant domain preserving discretization-independent schemes and convex limiting for hyperbolic systems.
+  [DOI: 10.1016/j.cma.2018.11.036](https://doi.org/10.1016/j.cma.2018.11.036)
+"""
+@inline function entropy_guermond_etal(u, equations::CompressibleEulerEquations3D)
+    rho, rho_v1, rho_v2, rho_v3, rho_e = u
+
+    # Modified specific entropy from Guermond et al. (2019)
+    s = (rho_e - 0.5f0 * (rho_v1^2 + rho_v2^2 + rho_v3^2) / rho) *
+        (1 / rho)^equations.gamma
+
+    return s
+end
+
+# Transformation from conservative variables u to d(s)/d(u)
+@inline function gradient_conservative(::typeof(entropy_guermond_etal),
+                                       u, equations::CompressibleEulerEquations3D)
+    return cons2entropy_guermond_etal(u, equations)
+end
+
 # Default entropy is the mathematical entropy
 @inline function entropy(cons, equations::CompressibleEulerEquations3D)
     return entropy_math(cons, equations)

src/solvers/dgsem_p4est/subcell_limiters_3d.jl

@@ -9,6 +9,187 @@
 # IDP Limiting
 ###############################################################################
 
+###############################################################################
+# Calculation of local bounds using low-order FV solution
+
+@inline function calc_bounds_onesided!(var_minmax, min_or_max, variable,
+                                       u::AbstractArray{<:Any, 5}, 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 k in eachnode(dg), j in eachnode(dg), i in eachnode(dg)
+            if min_or_max === max
+                var_minmax[i, j, k, element] = typemin(eltype(var_minmax))
+            else
+                var_minmax[i, j, k, element] = typemax(eltype(var_minmax))
+            end
+        end
+
+        # Calculate bounds at Gauss-Lobatto nodes
+        for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg)
+            var = variable(get_node_vars(u, equations, dg, i, j, k, element), equations)
+            var_minmax[i, j, k, element] = min_or_max(var_minmax[i, j, k, element], var)
+
+            if i > 1
+                var_minmax[i - 1, j, k, element] = min_or_max(var_minmax[i - 1, j, k,
+                                                                         element], var)
+            end
+            if i < nnodes(dg)
+                var_minmax[i + 1, j, k, element] = min_or_max(var_minmax[i + 1, j, k,
+                                                                         element], var)
+            end
+            if j > 1
+                var_minmax[i, j - 1, k, element] = min_or_max(var_minmax[i, j - 1, k,
+                                                                         element], var)
+            end
+            if j < nnodes(dg)
+                var_minmax[i, j + 1, k, element] = min_or_max(var_minmax[i, j + 1, k,
+                                                                         element], var)
+            end
+            if k > 1
+                var_minmax[i, j, k - 1, element] = min_or_max(var_minmax[i, j, k - 1,
+                                                                         element], var)
+            end
+            if k < nnodes(dg)
+                var_minmax[i, j, k + 1, element] = min_or_max(var_minmax[i, j, k + 1,
+                                                                         element], var)
+            end
+        end
+    end
+
+    # Values at element boundary
+    calc_bounds_onesided_interface!(var_minmax, min_or_max, variable, u, t, semi, mesh)
+
+    return nothing
+end
+
+function calc_bounds_onesided_interface!(var_minmax, minmax, variable, u, t, semi,
+                                         mesh::P4estMesh{3})
+    _, equations, dg, cache = mesh_equations_solver_cache(semi)
+    (; boundary_conditions) = semi
+
+    (; neighbor_ids, node_indices) = cache.interfaces
+    index_range = eachnode(dg)
+
+    # Calc bounds at interfaces and periodic boundaries
+    for interface in eachinterface(dg, cache)
+        # Get element and side index information on the primary element
+        primary_element = neighbor_ids[1, interface]
+        primary_indices = node_indices[1, interface]
+
+        # Get element and side index information on the secondary element
+        secondary_element = neighbor_ids[2, interface]
+        secondary_indices = node_indices[2, interface]
+
+        # Create the local i,j,k indexing
+        i_primary_start, i_primary_step_i, i_primary_step_j = index_to_start_step_3d(primary_indices[1],
+                                                                                     index_range)
+        j_primary_start, j_primary_step_i, j_primary_step_j = index_to_start_step_3d(primary_indices[2],
+                                                                                     index_range)
+        k_primary_start, k_primary_step_i, k_primary_step_j = index_to_start_step_3d(primary_indices[3],
+                                                                                     index_range)
+
+        i_primary = i_primary_start
+        j_primary = j_primary_start
+        k_primary = k_primary_start
+
+        i_secondary_start, i_secondary_step_i, i_secondary_step_j = index_to_start_step_3d(secondary_indices[1],
+                                                                                           index_range)
+        j_secondary_start, j_secondary_step_i, j_secondary_step_j = index_to_start_step_3d(secondary_indices[2],
+                                                                                           index_range)
+        k_secondary_start, k_secondary_step_i, k_secondary_step_j = index_to_start_step_3d(secondary_indices[3],
+                                                                                           index_range)
+
+        i_secondary = i_secondary_start
+        j_secondary = j_secondary_start
+        k_secondary = k_secondary_start
+
+        for j in eachnode(dg)
+            for i in eachnode(dg)
+                var_primary = variable(get_node_vars(u, equations, dg, i_primary,
+                                                     j_primary, k_primary,
+                                                     primary_element), equations)
+                var_secondary = variable(get_node_vars(u, equations, dg, i_secondary,
+                                                       j_secondary, k_secondary,
+                                                       secondary_element),
+                                         equations)
+
+                var_minmax[i_primary, j_primary, k_primary, primary_element] = minmax(var_minmax[i_primary,
+                                                                                                 j_primary,
+                                                                                                 k_primary,
+                                                                                                 primary_element],
+                                                                                      var_secondary)
+                var_minmax[i_secondary, j_secondary, k_secondary, secondary_element] = minmax(var_minmax[i_secondary,
+                                                                                                         j_secondary,
+                                                                                                         k_secondary,
+                                                                                                         secondary_element],
+                                                                                              var_primary)
+
+                # Increment the primary element indices
+                i_primary += i_primary_step_i
+                j_primary += j_primary_step_i
+                k_primary += k_primary_step_i
+                # Increment the secondary element surface indices
+                i_secondary += i_secondary_step_i
+                j_secondary += j_secondary_step_i
+                k_secondary += k_secondary_step_i
+            end
+            # Increment the primary element indices
+            i_primary += i_primary_step_j
+            j_primary += j_primary_step_j
+            k_primary += k_primary_step_j
+            # Increment the secondary element surface indices
+            i_secondary += i_secondary_step_j
+            j_secondary += j_secondary_step_j
+            k_secondary += k_secondary_step_j
+        end
+    end
+
+    # Calc bounds at physical boundaries
+    calc_bounds_onesided_boundary!(var_minmax, minmax, variable, u, t,
+                                   boundary_conditions,
+                                   mesh, equations, dg, cache)
+
+    return nothing
+end
+
+@inline function calc_bounds_onesided_boundary!(var_minmax, minmax, variable, u, t,
+                                                boundary_conditions::BoundaryConditionPeriodic,
+                                                mesh::P4estMesh{3},
+                                                equations, dg, cache)
+    return nothing
+end
+
+##############################################################################
+# Local one-sided limiting of nonlinear variables
+
+@inline function idp_local_onesided!(alpha, limiter, u::AbstractArray{<:Real, 5},
+                                     t, dt, semi,
+                                     variable, min_or_max)
+    mesh, equations, dg, cache = mesh_equations_solver_cache(semi)
+    (; variable_bounds) = limiter.cache.subcell_limiter_coefficients
+    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)
+        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)
+            u_local = get_node_vars(u, equations, dg, i, j, k, element)
+            newton_loops_alpha!(alpha, var_minmax[i, j, k, element],
+                                u_local, i, j, k, 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
+
+    return nothing
+end
+
 ###############################################################################
 # Global positivity limiting of conservative variables
 

test/test_p4est_3d.jl

@@ -371,21 +371,50 @@ end
     @test_allocations(Trixi.rhs!, semi, sol, 1000)
 end
 
-@trixi_testset "elixir_euler_sedov_sc_subcell.jl" begin
+@trixi_testset "elixir_euler_sedov_sc_subcell.jl (positivity bounds)" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_sedov_sc_subcell.jl"),
                         l2=[
-                            0.1942700455652903,
-                            0.07557644365785855,
-                            0.07557644365785836,
-                            0.07557644365785698,
-                            0.3713893635249306
+                            0.19427117014566905,
+                            0.07557679851688888,
+                            0.07557679851688896,
+                            0.07557679851688917,
+                            0.3713893373335463
                         ],
                         linf=[
-                            2.7542157588958798,
-                            1.8885700263691245,
-                            1.888570026369125,
-                            1.8885700263691252,
-                            4.9712792944452096
+                            2.754292081408987,
+                            1.888627709320322,
+                            1.8886277093203232,
+                            1.888627709320322,
+                            4.971280431903264
+                        ],
+                        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
+    # integrator which are not *recorded* for the methods from
+    # OrdinaryDiffEq.jl
+    # Corresponding issue: https://github.com/trixi-framework/Trixi.jl/issues/1877
+    @test_allocations(Trixi.rhs!, semi, sol, 15_000)
+end
+
+@trixi_testset "elixir_euler_sedov_sc_subcell.jl (local bounds)" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_sedov_sc_subcell.jl"),
+                        local_onesided_variables_nonlinear=[(entropy_guermond_etal,
+                                                             min)],
+                        max_iterations_newton=40,
+                        l2=[
+                            0.19153085678321066,
+                            0.07411109384422779,
+                            0.07411109384410808,
+                            0.07411109384406232,
+                            0.36714268468314665
+                        ],
+                        linf=[
+                            1.4037775549639524,
+                            1.339590863739464,
+                            1.3395908637591605,
+                            1.3395908637371077,
+                            4.824252924073932
                         ],
                         tspan=(0.0, 0.3))
     # Ensure that we do not have excessive memory allocations