Add one-sided local subcell IDP limiting for non-linear variables (tr…

…ixi-framework#1792)

* Add positivity limiting of non-linear variables

* Revise derivative function call; Add default derivative version

* Adapt test to actually test pos limiter for nonlinear variables

* Add unit test to test default implementation of variable_derivative

* Clean up comments and code

* Rename Newton-bisection variables

* Implement suggestions

* Relocate functions

* Add entropy limiters

* Update test errors after adding entropy limiting

* Fix bug in entropy limiting

* Adapt estimated errors after bug fix

* Remove doubled code

* Rename function

* 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

* Add comment to NEWS.md

* Remove whitespaces

* Implement suggestions

* Replace `foreach` with `for`

* Fix tests

* Use new bounds check reduction for one sided limiter

* Adapt tests after last merge of main

* Rename `entropy_spec` to `entropy_guermond_etal`

* Remove not-needed `tuple`

* Adapt nameing changes to tutorial

---------

Co-authored-by: Michael Schlottke-Lakemper <[email protected]>

NEWS.md

@@ -12,6 +12,7 @@ for human readability.
 - Implementation of 1D Linearized Euler Equations.
 - New analysis callback for 2D `P4estMesh` to compute integrated quantities along a boundary surface, e.g., pressure lift and drag coefficients.
 - Optional tuple parameter for `GlmSpeedCallback` called `semi_indices` to specify for which semidiscretization of a `SemidiscretizationCoupled` we need to update the GLM speed.
+- Subcell local one-sided limiting support for nonlinear variables in 2D for `TreeMesh`.
 
 ## Changes when updating to v0.7 from v0.6.x
 

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.

examples/tree_2d_dgsem/elixir_euler_blast_wave_sc_subcell_nonperiodic.jl

@@ -41,7 +41,9 @@ surface_flux = flux_lax_friedrichs
 volume_flux = flux_ranocha
 basis = LobattoLegendreBasis(3)
 limiter_idp = SubcellLimiterIDP(equations, basis;
-                                local_minmax_variables_cons = ["rho"])
+                                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)

examples/tree_2d_dgsem/elixir_euler_sedov_blast_wave_sc_subcell.jl

@@ -42,7 +42,9 @@ surface_flux = flux_lax_friedrichs
 volume_flux = flux_chandrashekar
 basis = LobattoLegendreBasis(3)
 limiter_idp = SubcellLimiterIDP(equations, basis;
-                                local_minmax_variables_cons = ["rho"])
+                                local_twosided_variables_cons = ["rho"],
+                                local_onesided_variables_nonlinear = [(Trixi.entropy_guermond_etal,
+                                                                       min)])
 volume_integral = VolumeIntegralSubcellLimiting(limiter_idp;
                                                 volume_flux_dg = volume_flux,
                                                 volume_flux_fv = surface_flux)

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)

src/callbacks_stage/subcell_bounds_check.jl

@@ -77,22 +77,27 @@ function init_callback(callback::BoundsCheckCallback, semi, limiter::SubcellLimi
         return nothing
     end
 
-    (; local_minmax, positivity) = 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 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")
@@ -120,15 +125,15 @@ end
 
 @inline function finalize_callback(callback::BoundsCheckCallback, semi,
                                    limiter::SubcellLimiterIDP)
-    (; local_minmax, positivity) = 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: ",
@@ -137,9 +142,19 @@ end
                     idp_bounds_delta_global[Symbol(v_string, "_max")])
         end
     end
+    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: ",

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_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
 
@@ -20,8 +20,8 @@
     # `@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_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")
@@ -45,9 +45,28 @@
             idp_bounds_delta_local[key_max] = deviation_max
         end
     end
+    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 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
+            end
+            idp_bounds_delta_local[key] = 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")
@@ -86,16 +105,23 @@
         # 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")],
                           ", ", idp_bounds_delta_local[Symbol(v_string, "_max")])
                 end
             end
+            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, ", ", idp_bounds_delta_local[Symbol(string(v), "_min")])

src/equations/compressible_euler_2d.jl

@@ -1886,6 +1886,27 @@ end
     return SVector(w1, w2, w3, w4)
 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::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
+
+    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)
@@ -1991,6 +2012,29 @@ end
     return S
 end
 
+# Transformation from conservative variables u to d(s)/d(u)
+@inline function gradient_conservative(::typeof(entropy_math),
+                                       u, equations::CompressibleEulerEquations2D)
+    return cons2entropy(u, equations)
+end
+
+# 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)
+    s = (rho_e - 0.5 * (rho_v1^2 + rho_v2^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::CompressibleEulerEquations2D)
+    return cons2entropy_guermond_etal(u, equations)
+end
+
 # Default entropy is the mathematical entropy
 @inline function entropy(cons, equations::CompressibleEulerEquations2D)
     entropy_math(cons, equations)

src/equations/ideal_glm_mhd_2d.jl

@@ -295,8 +295,8 @@ of local and symmetric parts. It is equivalent to the non-conservative flux of B
 et al. (`flux_nonconservative_powell`) for conforming meshes but it yields different
 results on non-conforming meshes(!).
 
-The two other flux functions with the same name return either the local 
-or symmetric portion of the non-conservative flux based on the type of the 
+The two other flux functions with the same name return either the local
+or symmetric portion of the non-conservative flux based on the type of the
 nonconservative_type argument, employing multiple dispatch. They are used to
 compute the subcell fluxes in dg_2d_subcell_limiters.jl.