Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Support nonlinear positivity limiting for arbitrary variables #115

Merged
merged 8 commits into from
Oct 17, 2023
Merged

Support nonlinear positivity limiting for arbitrary variables #115

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
a4a1fe8
Allow arbitrary nonlinear positivity limiting
bennibolm Oct 12, 2023
2eac4cf
Implement changes for limiting of entropys
bennibolm Oct 12, 2023
651b716
Relocate checks and goal functions
bennibolm Oct 12, 2023
70b4fd8
Rename functions
bennibolm Oct 12, 2023
2be3ba8
Merge branch 'subcell-limiting' into subcell-limiting-general-nonline…
bennibolm Oct 13, 2023
463de3a
Add @inline
bennibolm Oct 13, 2023
e4d3a9d
Add @inline to checks and goal functions; Update test required
bennibolm Oct 13, 2023
d04afad
Update tests
bennibolm Oct 13, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
7 changes: 6 additions & 1 deletion src/equations/compressible_euler_2d.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand All @@ -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)
Expand All @@ -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
Expand All @@ -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
Expand Down
111 changes: 54 additions & 57 deletions src/solvers/dgsem_tree/subcell_limiters_2d.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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, antidiffusive_flux2 = cache.antidiffusive_fluxes
if mesh isa TreeMesh
Expand All @@ -589,38 +564,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, 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, 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, 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, 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]
Expand All @@ -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
Expand All @@ -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
Expand All @@ -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
Expand All @@ -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
Expand All @@ -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)
Expand Down
6 changes: 3 additions & 3 deletions test/test_structured_2d.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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],
amrueda marked this conversation as resolved.
Show resolved Hide resolved
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

Expand Down Expand Up @@ -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
Expand Down