Skip to content

Commit

Permalink
Implement changes for limiting of entropys
Browse files Browse the repository at this point in the history
  • Loading branch information
@bennibolm
bennibolm committed Oct 12, 2023
1 parent a4a1fe8 commit 2eac4cf
Show file tree
Hide file tree
Showing 2 changed files with 31 additions and 34 deletions.
2 changes: 2 additions & 0 deletions src/equations/compressible_euler_2d.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1544,6 +1544,7 @@ end

return SVector(w1, w2, w3, w4)
end
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)
Expand All @@ -1570,6 +1571,7 @@ end

return SVector(w1, w2, w3, w4)
end
entropy_spec(u, equations, derivative::True) = cons2entropy_spec(u, equations)

# Transformation from conservative variables u to d(p)/d(u)
@inline function pressure(u, equations::CompressibleEulerEquations2D, derivative::True)
Expand Down
63 changes: 29 additions & 34 deletions src/solvers/dgsem_tree/subcell_limiters_2d.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -408,20 +408,16 @@ 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)
function initial_check_entropy_spec(bound, goal, newton_abstol)
goal <= max(newton_abstol, abs(bound) * newton_abstol)
end

Expand All @@ -438,20 +434,16 @@ 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)
function initial_check_entropy_math(bound, goal, newton_abstol)
goal >= -max(newton_abstol, abs(bound) * newton_abstol)
end

Expand Down Expand Up @@ -554,7 +546,9 @@ end

# Perform Newton's bisection method to find new alpha
newton_loops_alpha!(alpha, var_min[i, j, element], u_local, i, j, element,
variable, dt, mesh, equations, dg, cache, limiter)
variable, initial_check_nonnegative,
final_check_nonnegative,
dt, mesh, equations, dg, cache, limiter)
end
end

Expand All @@ -563,16 +557,17 @@ end

goal_function(variable, bound, u, equations) = bound - variable(u, equations)
function dgoal_function(variable, u, dt, antidiffusive_flux, equations)
-dot(variable(u, equations, true), dt * antidiffusive_flux)
-dot(variable(u, equations, True()), dt * antidiffusive_flux)
end

initialCheck(bound, goal, newton_abstol) = goal <= 0
function finalCheck(bound, goal, newton_abstol)
initial_check_nonnegative(bound, goal, newton_abstol) = goal <= 0
function final_check_nonnegative(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,
variable, 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, antidiffusive_flux2 = cache.antidiffusive_fluxes
if mesh isa TreeMesh
Expand All @@ -587,37 +582,37 @@ end
antidiffusive_flux = gamma_constant_newton * inverse_jacobian * inverse_weights[i] *
get_node_vars(antidiffusive_flux1, equations, dg, i, j,
element)
newton_loop!(alpha, bound, u, i, j, element, variable, 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, equations, dg, i + 1, j,
element)
newton_loop!(alpha, bound, u, i, j, element, variable, 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, equations, dg, i, j,
element)
newton_loop!(alpha, bound, u, i, j, element, variable, 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, equations, dg, i, j + 1,
element)
newton_loop!(alpha, bound, u, i, j, element, variable, 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, variable,
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]
Expand All @@ -631,7 +626,7 @@ end
if isValidState(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
Expand Down Expand Up @@ -666,7 +661,7 @@ end

# Check new beta for condition and update bounds
as = goal_function(variable, bound, u_curr, equations)
if initialCheck(bound, as, newton_abstol)
if initial_check(bound, as, newton_abstol)
# New beta fulfills condition
beta_L = beta
else
Expand All @@ -693,7 +688,7 @@ end
end

# Check absolute tolerance
if finalCheck(bound, as, newton_abstol)
if final_check(bound, as, newton_abstol)
break
end
end
Expand All @@ -708,7 +703,7 @@ end
return nothing
end

function standard_finalCheck(bound, goal, newton_abstol)
function final_check_standard(bound, goal, newton_abstol)
abs(goal) < max(newton_abstol, abs(bound) * newton_abstol)
end

Expand Down

0 comments on commit 2eac4cf

Please sign in to comment.