Merge branch 'subcell-limiting' into subcell_positivity_nonconservative

src/equations/compressible_euler_2d.jl

@@ -1544,6 +1544,7 @@ end
 
     return SVector(w1, w2, w3, w4)
 end
+@inline entropy_math(u, equations, derivative::True) = cons2entropy(u, equations)
 
 # 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)
@@ -1570,6 +1571,7 @@ end
 
     return SVector(w1, w2, w3, w4)
 end
+@inline entropy_spec(u, equations, derivative::True) = cons2entropy_spec(u, equations)
 
 # Transformation from conservative variables u to d(p)/d(u)
 @inline function dpdu(u, equations::CompressibleEulerEquations2D)
@@ -1581,6 +1583,9 @@ end
 
     return (equations.gamma - 1.0) * SVector(0.5 * v_square, -v1, -v2, 1.0)
 end
+@inline function pressure(u, equations::CompressibleEulerEquations2D, derivative::True)
+    return dpdu(u, equations)
+end
 
 @inline function entropy2cons(w, equations::CompressibleEulerEquations2D)
     # See Hughes, Franca, Mallet (1986) A new finite element formulation for CFD
@@ -1606,7 +1611,7 @@ end
     return SVector(rho, rho_v1, rho_v2, rho_e)
 end
 
-@inline function isValidState(cons, equations::CompressibleEulerEquations2D)
+@inline function is_valid_state(cons, equations::CompressibleEulerEquations2D)
     p = pressure(cons, equations)
     if cons[1] <= 0.0 || p <= 0.0
         return false

src/solvers/dgsem_tree/subcell_limiters_2d.jl

@@ -408,23 +408,15 @@ end
         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,
-                                specEntropy_goal, specEntropy_dGoal_dbeta,
-                                specEntropy_initialCheck, standard_finalCheck,
+                                entropy_spec, initial_check_entropy_spec,
+                                final_check_standard,
                                 dt, mesh, equations, dg, cache, limiter)
         end
     end
 
     return nothing
 end
 
-specEntropy_goal(bound, u, equations) = bound - entropy_spec(u, equations)
-function specEntropy_dGoal_dbeta(u, dt, antidiffusive_flux, equations)
-    -dot(cons2entropy_spec(u, equations), dt * antidiffusive_flux)
-end
-function specEntropy_initialCheck(bound, goal, newton_abstol)
-    goal <= max(newton_abstol, abs(bound) * newton_abstol)
-end
-
 @inline function idp_math_entropy!(alpha, limiter, u, t, dt, semi, elements)
     mesh, equations, dg, cache = mesh_equations_solver_cache(semi)
     (; variable_bounds) = limiter.cache.subcell_limiter_coefficients
@@ -438,23 +430,15 @@ end
         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,
-                                mathEntropy_goal, mathEntropy_dGoal_dbeta,
-                                mathEntropy_initialCheck, standard_finalCheck,
+                                entropy_math, initial_check_entropy_math,
+                                final_check_standard,
                                 dt, mesh, equations, dg, cache, limiter)
         end
     end
 
     return nothing
 end
 
-mathEntropy_goal(bound, u, equations) = bound - entropy_math(u, equations)
-function mathEntropy_dGoal_dbeta(u, dt, antidiffusive_flux, equations)
-    -dot(cons2entropy(u, equations), dt * antidiffusive_flux)
-end
-function mathEntropy_initialCheck(bound, goal, newton_abstol)
-    goal >= -max(newton_abstol, abs(bound) * newton_abstol)
-end
-
 @inline function idp_positivity!(alpha, limiter, u, dt, semi, elements)
     # Conservative variables
     for variable in limiter.positivity_variables_cons
@@ -554,27 +538,18 @@ end
 
             # Perform Newton's bisection method to find new alpha
             newton_loops_alpha!(alpha, var_min[i, j, element], u_local, i, j, element,
-                                pressure_goal, pressure_dgoal_dbeta,
-                                pressure_initialCheck, pressure_finalCheck,
+                                variable, initial_check_nonnegative,
+                                final_check_nonnegative,
                                 dt, mesh, equations, dg, cache, limiter)
         end
     end
 
     return nothing
 end
 
-pressure_goal(bound, u, equations) = bound - pressure(u, equations)
-function pressure_dgoal_dbeta(u, dt, antidiffusive_flux, equations)
-    -dot(dpdu(u, equations), dt * antidiffusive_flux)
-end
-pressure_initialCheck(bound, goal, newton_abstol) = goal <= 0
-function pressure_finalCheck(bound, goal, newton_abstol)
-    (goal <= eps()) && (goal > -max(newton_abstol, abs(bound) * newton_abstol))
-end
-
-@inline function newton_loops_alpha!(alpha, bound, u, i, j, element,
-                                     goal_fct, dgoal_fct, initialCheck, finalCheck,
-                                     dt, mesh, equations, dg, cache, limiter)
+@inline function newton_loops_alpha!(alpha, bound, u, i, j, element, variable,
+                                     initial_check, final_check, dt, mesh, equations,
+                                     dg, cache, limiter)
     @unpack inverse_weights = dg.basis
     @unpack antidiffusive_flux1_L, antidiffusive_flux2_L, antidiffusive_flux1_R, antidiffusive_flux2_R = cache.antidiffusive_fluxes
     if mesh isa TreeMesh
@@ -589,38 +564,37 @@ 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, goal_fct, dgoal_fct, initialCheck,
-                 finalCheck, equations, dt, limiter, antidiffusive_flux)
+    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, goal_fct, dgoal_fct, initialCheck,
-                 finalCheck, equations, dt, limiter, antidiffusive_flux)
+    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, goal_fct, dgoal_fct, initialCheck,
-                 finalCheck, equations, dt, limiter, antidiffusive_flux)
+    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, goal_fct, dgoal_fct, initialCheck,
-                 finalCheck, equations, dt, limiter, antidiffusive_flux)
+    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,
-                              goal_fct, dgoal_fct, initialCheck, finalCheck,
-                              equations, dt, limiter, antidiffusive_flux)
+@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]
@@ -631,19 +605,20 @@ end
     u_curr = u + beta * dt * antidiffusive_flux
 
     # If state is valid, perform initial check and return if correction is not needed
-    if isValidState(u_curr, equations)
-        as = goal_fct(bound, u_curr, equations)
+    if is_valid_state(u_curr, equations)
+        as = goal_function(variable, bound, u_curr, equations)
 
-        initialCheck(bound, as, newton_abstol) && return nothing
+        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 isValidState(u_curr, equations)
-            dSdbeta = dgoal_fct(u_curr, dt, antidiffusive_flux, equations)
+        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
@@ -661,14 +636,14 @@ end
             u_curr = u + beta * dt * antidiffusive_flux
 
             # If the state is invalid, finish bisection step without checking tolerance and iterate further
-            if !isValidState(u_curr, equations)
+            if !is_valid_state(u_curr, equations)
                 beta_R = beta
                 continue
             end
 
             # Check new beta for condition and update bounds
-            as = goal_fct(bound, u_curr, equations)
-            if initialCheck(bound, as, newton_abstol)
+            as = goal_function(variable, bound, u_curr, equations)
+            if initial_check(bound, as, newton_abstol)
                 # New beta fulfills condition
                 beta_L = beta
             else
@@ -680,13 +655,13 @@ end
             u_curr = u + beta * dt * antidiffusive_flux
 
             # If the state is invalid, redefine right bound without checking tolerance and iterate further
-            if !isValidState(u_curr, equations)
+            if !is_valid_state(u_curr, equations)
                 beta_R = beta
                 continue
             end
 
             # Evaluate goal function
-            as = goal_fct(bound, u_curr, equations)
+            as = goal_function(variable, bound, u_curr, equations)
         end
 
         # Check relative tolerance
@@ -695,7 +670,7 @@ end
         end
 
         # Check absolute tolerance
-        if finalCheck(bound, as, newton_abstol)
+        if final_check(bound, as, newton_abstol)
             break
         end
     end
@@ -710,10 +685,32 @@ end
     return nothing
 end
 
-function standard_finalCheck(bound, goal, newton_abstol)
+# 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
+@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 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
+
 # this method is used when the limiter is constructed as for shock-capturing volume integrals
 function create_cache(limiter::Type{SubcellLimiterMCL}, equations::AbstractEquations{2},
                       basis::LobattoLegendreBasis, PressurePositivityLimiterKuzmin)

test/test_structured_2d.jl

@@ -187,8 +187,8 @@ isdir(outdir) && rm(outdir, recursive=true)
 
   @trixi_testset "elixir_euler_source_terms_sc_subcell.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms_sc_subcell.jl"),
-      l2   = [0.008160127272557726, 0.008658253869683077, 0.009351900401871649, 0.02775701488343099],
-      linf = [0.027225608222781528, 0.0407340321806311, 0.0381940733564341, 0.08080650914262844],
+      l2   = [0.00816013114351954, 0.008658251709937477, 0.009351905651482216, 0.027757012781694318],
+      linf = [0.027225615981281148, 0.040734036539016305, 0.0381940733564341, 0.08080650914262844],
       tspan = (0.0, 0.5))
   end
 
@@ -248,7 +248,7 @@ isdir(outdir) && rm(outdir, recursive=true)
   @trixi_testset "elixir_euler_shock_upstream_sc_subcell.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_shock_upstream_sc_subcell.jl"),
       l2   = [1.2351468819080416, 1.1269856120551724, 1.7239124305681928, 11.715260007491556],
-      linf = [5.385491808683259, 6.575446013701839, 10.065227889186632, 51.008985921289565],
+      linf = [5.385492532917423, 6.575446146030286, 10.0652310822613, 51.00901293102744],
       cells_per_dimension = (8, 12),
       tspan = (0.0, 0.5))
   end