Skip to content

Commit

Permalink
Add subcell positivity limiting of non-linear variables (#1738)
Browse files Browse the repository at this point in the history
* 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

* Implement suggestions

* Change error message for negative value with low-order method

* Add changes from main to new limiter

* Update NEWS.md

* Rename is_valid_state and gradient_u

---------

Co-authored-by: Michael Schlottke-Lakemper <[email protected]>
  • Loading branch information
bennibolm and sloede authored Jan 31, 2024
1 parent c7cee98 commit dfd632e
Show file tree
Hide file tree
Showing 13 changed files with 496 additions and 35 deletions.
1 change: 1 addition & 0 deletions NEWS.md
Original file line number Diff line number Diff line change
Expand Up @@ -11,6 +11,7 @@ for human readability.
- `flux_hllc` on non-cartesian meshes for `CompressibleEulerEquations{2,3}D`
- Different boundary conditions for quad/hex meshes in Abaqus format, even if not generated by HOHQMesh,
can now be digested by Trixi in 2D and 3D.
- Subcell (positivity) limiting support for nonlinear variables in 2D for `TreeMesh`

## Changes when updating to v0.6 from v0.5.x

Expand Down
Original file line number Diff line number Diff line change
@@ -0,0 +1,91 @@

using OrdinaryDiffEq
using Trixi

###############################################################################
# semidiscretization of the compressible Euler equations
gamma = 1.4
equations = CompressibleEulerEquations2D(gamma)

"""
initial_condition_kelvin_helmholtz_instability(x, t, equations::CompressibleEulerEquations2D)
A version of the classical Kelvin-Helmholtz instability based on
- Andrés M. Rueda-Ramírez, Gregor J. Gassner (2021)
A Subcell Finite Volume Positivity-Preserving Limiter for DGSEM Discretizations
of the Euler Equations
[arXiv: 2102.06017](https://arxiv.org/abs/2102.06017)
"""
function initial_condition_kelvin_helmholtz_instability(x, t,
equations::CompressibleEulerEquations2D)
# change discontinuity to tanh
# typical resolution 128^2, 256^2
# domain size is [-1,+1]^2
slope = 15
amplitude = 0.02
B = tanh(slope * x[2] + 7.5) - tanh(slope * x[2] - 7.5)
rho = 0.5 + 0.75 * B
v1 = 0.5 * (B - 1)
v2 = 0.1 * sin(2 * pi * x[1])
p = 1.0
return prim2cons(SVector(rho, v1, v2, p), equations)
end
initial_condition = initial_condition_kelvin_helmholtz_instability

surface_flux = flux_lax_friedrichs
volume_flux = flux_ranocha
polydeg = 3
basis = LobattoLegendreBasis(polydeg)

limiter_idp = SubcellLimiterIDP(equations, basis;
positivity_variables_cons = ["rho"],
positivity_variables_nonlinear = [pressure])
volume_integral = VolumeIntegralSubcellLimiting(limiter_idp;
volume_flux_dg = volume_flux,
volume_flux_fv = surface_flux)
solver = DGSEM(basis, surface_flux, volume_integral)

coordinates_min = (-1.0, -1.0)
coordinates_max = (1.0, 1.0)
mesh = TreeMesh(coordinates_min, coordinates_max,
initial_refinement_level = 5,
n_cells_max = 100_000)
semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)

###############################################################################
# ODE solvers, callbacks etc.

tspan = (0.0, 3.7)
ode = semidiscretize(semi, tspan)

summary_callback = SummaryCallback()

analysis_interval = 1000
analysis_callback = AnalysisCallback(semi, interval = analysis_interval)

alive_callback = AliveCallback(analysis_interval = analysis_interval)

save_solution = SaveSolutionCallback(interval = 100,
save_initial_solution = true,
save_final_solution = true,
solution_variables = cons2prim)

save_restart = SaveRestartCallback(interval = 1000,
save_final_restart = true)

stepsize_callback = StepsizeCallback(cfl = 0.7)

callbacks = CallbackSet(summary_callback,
analysis_callback, alive_callback,
stepsize_callback,
save_restart, save_solution)

###############################################################################
# run the simulation

stage_callbacks = (SubcellLimiterIDPCorrection(), BoundsCheckCallback(save_errors = false))

sol = Trixi.solve(ode, Trixi.SimpleSSPRK33(stage_callbacks = stage_callbacks);
dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
callback = callbacks);
summary_callback() # print the timer summary
7 changes: 4 additions & 3 deletions examples/tree_2d_dgsem/elixir_mhd_shockcapturing_subcell.jl
Original file line number Diff line number Diff line change
Expand Up @@ -22,7 +22,7 @@ function initial_condition_blast_wave(x, t, equations::IdealGlmMhdEquations2D)
r = sqrt(x[1]^2 + x[2]^2)

pmax = 10.0
pmin = 1.0
pmin = 0.01
rhomax = 1.0
rhomin = 0.01
if r <= 0.09
Expand Down Expand Up @@ -52,7 +52,8 @@ basis = LobattoLegendreBasis(3)

limiter_idp = SubcellLimiterIDP(equations, basis;
positivity_variables_cons = ["rho"],
positivity_correction_factor = 0.5)
positivity_variables_nonlinear = [pressure],
positivity_correction_factor = 0.1)
volume_integral = VolumeIntegralSubcellLimiting(limiter_idp;
volume_flux_dg = volume_flux,
volume_flux_fv = surface_flux)
Expand Down Expand Up @@ -84,7 +85,7 @@ save_solution = SaveSolutionCallback(interval = 100,
save_final_solution = true,
solution_variables = cons2prim)

cfl = 0.5
cfl = 0.4
stepsize_callback = StepsizeCallback(cfl = cfl)

glm_speed_callback = GlmSpeedCallback(glm_scale = 0.5, cfl = cfl)
Expand Down
8 changes: 8 additions & 0 deletions src/callbacks_stage/subcell_bounds_check.jl
Original file line number Diff line number Diff line change
Expand Up @@ -97,6 +97,9 @@ function init_callback(callback::BoundsCheckCallback, semi, limiter::SubcellLimi
end
print(f, ", " * string(variables[v]) * "_min")
end
for variable in limiter.positivity_variables_nonlinear
print(f, ", " * string(variable) * "_min")
end
end
println(f)
end
Expand Down Expand Up @@ -142,6 +145,11 @@ end
println(string(variables[v]) * ":\n- positivity: ",
idp_bounds_delta_global[Symbol(string(v), "_min")])
end
for variable in limiter.positivity_variables_nonlinear
variable_string = string(variable)
println(variable_string * ":\n- positivity: ",
idp_bounds_delta_global[Symbol(variable_string, "_min")])
end
end
println(""^100 * "\n")

Expand Down
18 changes: 18 additions & 0 deletions src/callbacks_stage/subcell_bounds_check_2d.jl
Original file line number Diff line number Diff line change
Expand Up @@ -60,6 +60,20 @@
deviation_threaded[stride_size * Threads.threadid()] = deviation
end
end
for variable in limiter.positivity_variables_nonlinear
key = Symbol(string(variable), "_min")
deviation_threaded = idp_bounds_delta_local[key]
@threaded for element in eachelement(solver, cache)
deviation = deviation_threaded[stride_size * Threads.threadid()]
for j in eachnode(solver), i in eachnode(solver)
var = variable(get_node_vars(u, equations, solver, i, j, element),
equations)
deviation = max(deviation,
variable_bounds[key][i, j, element] - var)
end
deviation_threaded[stride_size * Threads.threadid()] = deviation
end
end
end

for (key, _) in idp_bounds_delta_local
Expand Down Expand Up @@ -92,6 +106,10 @@
print(f, ", ",
idp_bounds_delta_local[Symbol(string(v), "_min")][stride_size])
end
for variable in limiter.positivity_variables_nonlinear
print(f, ", ",
idp_bounds_delta_local[Symbol(string(variable), "_min")][stride_size])
end
end
println(f)
end
Expand Down
21 changes: 21 additions & 0 deletions src/equations/compressible_euler_2d.jl
Original file line number Diff line number Diff line change
Expand Up @@ -1632,6 +1632,18 @@ end
return p
end

# Transformation from conservative variables u to d(p)/d(u)
@inline function gradient_conservative(::typeof(pressure),
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

return (equations.gamma - 1.0) * SVector(0.5 * v_square, -v1, -v2, 1.0)
end

@inline function density_pressure(u, equations::CompressibleEulerEquations2D)
rho, rho_v1, rho_v2, rho_e = u
rho_times_p = (equations.gamma - 1) * (rho * rho_e - 0.5 * (rho_v1^2 + rho_v2^2))
Expand Down Expand Up @@ -1699,4 +1711,13 @@ end
@inline function energy_internal(cons, equations::CompressibleEulerEquations2D)
return energy_total(cons, equations) - energy_kinetic(cons, equations)
end

# State validation for Newton-bisection method of subcell IDP limiting
@inline function Base.isvalid(u, equations::CompressibleEulerEquations2D)
p = pressure(u, equations)
if u[1] <= 0.0 || p <= 0.0
return false
end
return true
end
end # @muladd
6 changes: 6 additions & 0 deletions src/equations/equations.jl
Original file line number Diff line number Diff line change
Expand Up @@ -376,6 +376,12 @@ of the correct length `nvariables(equations)`.
"""
function energy_internal end

# Default implementation of gradient for `variable`. Used for subcell limiting.
# Implementing a gradient function for a specific variable improves the performance.
@inline function gradient_conservative(variable, u, equations)
return ForwardDiff.gradient(x -> variable(x, equations), u)
end

####################################################################################################
# Include files with actual implementations for different systems of equations.

Expand Down
23 changes: 23 additions & 0 deletions src/equations/ideal_glm_mhd_2d.jl
Original file line number Diff line number Diff line change
Expand Up @@ -1118,6 +1118,20 @@ end
return p
end

# Transformation from conservative variables u to d(p)/d(u)
@inline function gradient_conservative(::typeof(pressure),
u, equations::IdealGlmMhdEquations2D)
rho, rho_v1, rho_v2, rho_v3, rho_e, B1, B2, B3, psi = u

v1 = rho_v1 / rho
v2 = rho_v2 / rho
v3 = rho_v3 / rho
v_square = v1^2 + v2^2 + v3^2

return (equations.gamma - 1.0) *
SVector(0.5 * v_square, -v1, -v2, -v3, 1.0, -B1, -B2, -B3, -psi)
end

@inline function density_pressure(u, equations::IdealGlmMhdEquations2D)
rho, rho_v1, rho_v2, rho_v3, rho_e, B1, B2, B3, psi = u
p = (equations.gamma - 1) * (rho_e - 0.5 * (rho_v1^2 + rho_v2^2 + rho_v3^2) / rho
Expand Down Expand Up @@ -1384,6 +1398,15 @@ end
cons[9]^2 / 2)
end

# State validation for Newton-bisection method of subcell IDP limiting
@inline function Base.isvalid(u, equations::IdealGlmMhdEquations2D)
p = pressure(u, equations)
if u[1] <= 0.0 || p <= 0.0
return false
end
return true
end

# Calculate the cross helicity (\vec{v}⋅\vec{B}) for a conservative state `cons'
@inline function cross_helicity(cons, ::IdealGlmMhdEquations2D)
return (cons[2] * cons[6] + cons[3] * cons[7] + cons[4] * cons[8]) / cons[1]
Expand Down
60 changes: 46 additions & 14 deletions src/solvers/dgsem_tree/subcell_limiters.jl
Original file line number Diff line number Diff line change
Expand Up @@ -16,18 +16,28 @@ end
SubcellLimiterIDP(equations::AbstractEquations, basis;
local_minmax_variables_cons = String[],
positivity_variables_cons = String[],
positivity_correction_factor = 0.1)
positivity_variables_nonlinear = [],
positivity_correction_factor = 0.1,
max_iterations_newton = 10,
newton_tolerances = (1.0e-12, 1.0e-14),
gamma_constant_newton = 2 * ndims(equations))
Subcell invariant domain preserving (IDP) limiting used with [`VolumeIntegralSubcellLimiting`](@ref)
including:
- Local maximum/minimum Zalesak-type limiting for conservative variables (`local_minmax_variables_cons`)
- Positivity limiting for conservative variables (`positivity_variables_cons`)
- Positivity limiting for conservative variables (`positivity_variables_cons`) and nonlinear variables
(`positivity_variables_nonlinear`)
Conservative variables to be limited are passed as a vector of strings, e.g. `local_minmax_variables_cons = ["rho"]`
and `positivity_variables_cons = ["rho"]`.
and `positivity_variables_cons = ["rho"]`. For nonlinear variables the specific functions are
passed in a vector, e.g. `positivity_variables_nonlinear = [pressure]`.
The bounds are calculated using the low-order FV solution. The positivity limiter uses
`positivity_correction_factor` such that `u^new >= positivity_correction_factor * u^FV`.
The limiting of nonlinear variables uses a Newton-bisection method with a maximum of
`max_iterations_newton` iterations, relative and absolute tolerances of `newton_tolerances`
and a provisional update constant `gamma_constant_newton` (`gamma_constant_newton>=2*d`,
where `d = #dimensions`). See equation (20) of Pazner (2020) and equation (30) of Rueda-Ramírez et al. (2022).
!!! note
This limiter and the correction callback [`SubcellLimiterIDPCorrection`](@ref) only work together.
Expand All @@ -45,22 +55,32 @@ The bounds are calculated using the low-order FV solution. The positivity limite
!!! warning "Experimental implementation"
This is an experimental feature and may change in future releases.
"""
struct SubcellLimiterIDP{RealT <: Real, Cache} <: AbstractSubcellLimiter
struct SubcellLimiterIDP{RealT <: Real, LimitingVariablesNonlinear, Cache} <:
AbstractSubcellLimiter
local_minmax::Bool
local_minmax_variables_cons::Vector{Int} # Local mininum/maximum principles for conservative variables
positivity::Bool
positivity_variables_cons::Vector{Int} # Positivity for conservative variables
positivity_variables_nonlinear::LimitingVariablesNonlinear # Positivity for nonlinear variables
positivity_correction_factor::RealT
cache::Cache
max_iterations_newton::Int
newton_tolerances::Tuple{RealT, RealT} # Relative and absolute tolerances for Newton's method
gamma_constant_newton::RealT # Constant for the subcell limiting of convex (nonlinear) constraints
end

# this method is used when the limiter is constructed as for shock-capturing volume integrals
function SubcellLimiterIDP(equations::AbstractEquations, basis;
local_minmax_variables_cons = String[],
positivity_variables_cons = String[],
positivity_correction_factor = 0.1)
positivity_variables_nonlinear = [],
positivity_correction_factor = 0.1,
max_iterations_newton = 10,
newton_tolerances = (1.0e-12, 1.0e-14),
gamma_constant_newton = 2 * ndims(equations))
local_minmax = (length(local_minmax_variables_cons) > 0)
positivity = (length(positivity_variables_cons) > 0)
positivity = (length(positivity_variables_cons) +
length(positivity_variables_nonlinear) > 0)

local_minmax_variables_cons_ = get_variable_index.(local_minmax_variables_cons,
equations)
Expand All @@ -80,13 +100,20 @@ function SubcellLimiterIDP(equations::AbstractEquations, basis;
bound_keys = (bound_keys..., Symbol(string(v), "_min"))
end
end
for variable in positivity_variables_nonlinear
bound_keys = (bound_keys..., Symbol(string(variable), "_min"))
end

cache = create_cache(SubcellLimiterIDP, equations, basis, bound_keys)

SubcellLimiterIDP{typeof(positivity_correction_factor),
typeof(positivity_variables_nonlinear),
typeof(cache)}(local_minmax, local_minmax_variables_cons_,
positivity, positivity_variables_cons_,
positivity_correction_factor, cache)
positivity_variables_nonlinear,
positivity_correction_factor, cache,
max_iterations_newton, newton_tolerances,
gamma_constant_newton)
end

function Base.show(io::IO, limiter::SubcellLimiterIDP)
Expand All @@ -97,10 +124,15 @@ function Base.show(io::IO, limiter::SubcellLimiterIDP)
if !(local_minmax || positivity)
print(io, "No limiter selected => pure DG method")
else
print(io, "limiter=(")
local_minmax && print(io, "min/max limiting, ")
positivity && print(io, "positivity")
print(io, "), ")
features = String[]
if local_minmax
push!(features, "local min/max")
end
if positivity
push!(features, "positivity")
end
join(io, features, ", ")
print(io, "Limiter=($features), ")
end
print(io, "Local bounds with FV solution")
print(io, ")")
Expand All @@ -120,15 +152,15 @@ function Base.show(io::IO, ::MIME"text/plain", limiter::SubcellLimiterIDP)
if local_minmax
setup = [
setup...,
"" => "local maximum/minimum bounds for conservative variables $(limiter.local_minmax_variables_cons)",
"" => "Local maximum/minimum limiting for conservative variables $(limiter.local_minmax_variables_cons)",
]
end
if positivity
string = "positivity for conservative variables $(limiter.positivity_variables_cons)"
string = "Positivity limiting for conservative variables $(limiter.positivity_variables_cons) and $(limiter.positivity_variables_nonlinear)"
setup = [setup..., "" => string]
setup = [
setup...,
"" => " positivity correction factor = $(limiter.positivity_correction_factor)",
"" => "- with positivity correction factor = $(limiter.positivity_correction_factor)",
]
end
setup = [
Expand Down
Loading

0 comments on commit dfd632e

Please sign in to comment.