diff --git a/src/equations/compressible_euler_2d.jl b/src/equations/compressible_euler_2d.jl index d95fbdbcfc9..8fd99bbfc6b 100644 --- a/src/equations/compressible_euler_2d.jl +++ b/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 diff --git a/src/solvers/dgsem_tree/subcell_limiters_2d.jl b/src/solvers/dgsem_tree/subcell_limiters_2d.jl index c10c6512b83..4290c327ad5 100644 --- a/src/solvers/dgsem_tree/subcell_limiters_2d.jl +++ b/src/solvers/dgsem_tree/subcell_limiters_2d.jl @@ -408,8 +408,8 @@ 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 @@ -417,14 +417,6 @@ 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,8 +430,8 @@ 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 @@ -447,14 +439,6 @@ 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,8 +538,8 @@ 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 @@ -563,18 +547,9 @@ 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,10 +605,10 @@ 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 @@ -642,8 +616,9 @@ end 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) diff --git a/test/test_structured_2d.jl b/test/test_structured_2d.jl index be1a1d3138b..f776de954bc 100644 --- a/test/test_structured_2d.jl +++ b/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