Skip to content

Add local limiting support for nonperiodic TreeMeshes #2878

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

Open
bennibolm wants to merge 7 commits into
base: main
Choose a base branch
from
Open

Add local limiting support for nonperiodic TreeMeshes #2878

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
621df7e
Allow local limiting for nonperiodic TreeMeshes
bennibolm Mar 23, 2026
3c85c9b
Merge branch 'main' into subcell-3d-nonperiodic-local-limiting
bennibolm Mar 30, 2026
ece9345
Introduce `calc_bounds_one/twosided_boundary` for Structured and Tree…
bennibolm Mar 30, 2026
a2056b3
Fix ambiguity
bennibolm Mar 31, 2026
7813812
Merge branch 'main' into subcell-3d-nonperiodic-local-limiting
DanielDoehring Mar 31, 2026
da6b7db
Merge branch 'main' into subcell-3d-nonperiodic-local-limiting
DanielDoehring Mar 31, 2026
e981a7b
Merge branch 'main' into subcell-3d-nonperiodic-local-limiting
bennibolm Apr 2, 2026
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
4 changes: 2 additions & 2 deletions examples/tree_3d_dgsem/elixir_euler_sedov_blast_wave_sc_subcell.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -59,11 +59,11 @@ coordinates_max = (1.0, 1.0, 1.0)
mesh = TreeMesh(coordinates_min, coordinates_max,
initial_refinement_level = 3,
n_cells_max = 100_000,
periodicity = true)
periodicity = false)

# create the semi discretization object
semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
boundary_conditions = boundary_condition_periodic)
boundary_conditions = BoundaryConditionDirichlet(initial_condition))

###############################################################################
# ODE solvers, callbacks etc.
Expand Down
5 changes: 2 additions & 3 deletions src/solvers/dgsem_p4est/subcell_limiters_2d.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -8,7 +8,6 @@
function calc_bounds_twosided_interface!(var_min, var_max, variable,
u, t, semi, mesh::P4estMesh{2}, equations)
_, _, dg, cache = mesh_equations_solver_cache(semi)
(; boundary_conditions) = semi

(; neighbor_ids, node_indices) = cache.interfaces
index_range = eachnode(dg)
Expand Down Expand Up @@ -69,6 +68,7 @@ function calc_bounds_twosided_interface!(var_min, var_max, variable,
end

# Calc bounds at physical boundaries
(; boundary_conditions) = semi
calc_bounds_twosided_boundary!(var_min, var_max, variable, u, t,
boundary_conditions,
mesh, equations, dg, cache)
Expand Down Expand Up @@ -136,11 +136,9 @@ end
function calc_bounds_onesided_interface!(var_minmax, minmax, variable, u, t, semi,
mesh::P4estMesh{2})
_, equations, dg, cache = mesh_equations_solver_cache(semi)
(; boundary_conditions) = semi

(; neighbor_ids, node_indices) = cache.interfaces
index_range = eachnode(dg)
index_end = last(index_range)

# Calc bounds at interfaces and periodic boundaries
for interface in eachinterface(dg, cache)
Expand Down Expand Up @@ -192,6 +190,7 @@ function calc_bounds_onesided_interface!(var_minmax, minmax, variable, u, t, sem
end

# Calc bounds at physical boundaries
(; boundary_conditions) = semi
calc_bounds_onesided_boundary!(var_minmax, minmax, variable, u, t,
boundary_conditions,
mesh, equations, dg, cache)
Expand Down
4 changes: 2 additions & 2 deletions src/solvers/dgsem_p4est/subcell_limiters_3d.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -8,7 +8,6 @@
function calc_bounds_twosided_interface!(var_min, var_max, variable,
u, t, semi, mesh::P4estMesh{3}, equations)
_, _, dg, cache = mesh_equations_solver_cache(semi)
(; boundary_conditions) = semi

(; neighbor_ids, node_indices) = cache.interfaces
index_range = eachnode(dg)
Expand Down Expand Up @@ -96,6 +95,7 @@ function calc_bounds_twosided_interface!(var_min, var_max, variable,
end

# Calc bounds at physical boundaries
(; boundary_conditions) = semi
calc_bounds_twosided_boundary!(var_min, var_max, variable, u, t,
boundary_conditions,
mesh, equations, dg, cache)
Expand Down Expand Up @@ -180,7 +180,6 @@ end
function calc_bounds_onesided_interface!(var_minmax, minmax, variable, u, t, semi,
mesh::P4estMesh{3})
_, equations, dg, cache = mesh_equations_solver_cache(semi)
(; boundary_conditions) = semi

(; neighbor_ids, node_indices) = cache.interfaces
index_range = eachnode(dg)
Expand Down Expand Up @@ -260,6 +259,7 @@ function calc_bounds_onesided_interface!(var_minmax, minmax, variable, u, t, sem
end

# Calc bounds at physical boundaries
(; boundary_conditions) = semi
calc_bounds_onesided_boundary!(var_minmax, minmax, variable, u, t,
boundary_conditions,
mesh, equations, dg, cache)
Expand Down
52 changes: 42 additions & 10 deletions src/solvers/dgsem_structured/subcell_limiters_2d.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -8,8 +8,6 @@
function calc_bounds_twosided_interface!(var_min, var_max, variable,
u, t, semi, mesh::StructuredMesh{2}, equations)
_, _, dg, cache = mesh_equations_solver_cache(semi)
(; boundary_conditions) = semi
(; contravariant_vectors) = cache.elements

# Calc bounds at interfaces and periodic boundaries
for element in eachelement(dg, cache)
Expand Down Expand Up @@ -48,9 +46,27 @@ function calc_bounds_twosided_interface!(var_min, var_max, variable,
end

# Calc bounds at physical boundaries
if isperiodic(mesh)
return nothing
end
(; boundary_conditions) = semi
calc_bounds_twosided_boundary!(var_min, var_max, variable, u, t,
boundary_conditions,
mesh, equations, dg, cache)

return nothing
end

@inline function calc_bounds_twosided_boundary!(var_min, var_max, variable, u, t,
boundary_conditions::BoundaryConditionPeriodic,
mesh::StructuredMesh{2},
equations, dg, cache)
return nothing
end

@inline function calc_bounds_twosided_boundary!(var_min, var_max, variable, u, t,
boundary_conditions,
mesh::StructuredMesh{2},
equations, dg, cache)
(; contravariant_vectors) = cache.elements

linear_indices = LinearIndices(size(mesh))
if !isperiodic(mesh, 1)
# - xi direction
Expand Down Expand Up @@ -135,8 +151,6 @@ end
function calc_bounds_onesided_interface!(var_minmax, minmax, variable, u, t, semi,
mesh::StructuredMesh{2})
_, equations, dg, cache = mesh_equations_solver_cache(semi)
(; boundary_conditions) = semi
(; contravariant_vectors) = cache.elements

# Calc bounds at interfaces and periodic boundaries
for element in eachelement(dg, cache)
Expand Down Expand Up @@ -172,9 +186,27 @@ function calc_bounds_onesided_interface!(var_minmax, minmax, variable, u, t, sem
end

# Calc bounds at physical boundaries
if isperiodic(mesh)
return nothing
end
(; boundary_conditions) = semi
calc_bounds_onesided_boundary!(var_minmax, minmax, variable, u, t,
boundary_conditions,
mesh, equations, dg, cache)

return nothing
end

@inline function calc_bounds_onesided_boundary!(var_minmax, minmax, variable, u, t,
boundary_conditions::BoundaryConditionPeriodic,
mesh::StructuredMesh{2},
equations, dg, cache)
return nothing
end

@inline function calc_bounds_onesided_boundary!(var_minmax, minmax, variable, u, t,
boundary_conditions,
mesh::StructuredMesh{2},
equations, dg, cache)
(; contravariant_vectors) = cache.elements

linear_indices = LinearIndices(size(mesh))
if !isperiodic(mesh, 1)
# - xi direction
Expand Down
42 changes: 40 additions & 2 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 @@ -57,7 +57,7 @@ end
u, t, semi, mesh::TreeMesh2D,
equations)
_, _, dg, cache = mesh_equations_solver_cache(semi)
(; boundary_conditions) = semi

# Calc bounds at interfaces and periodic boundaries
for interface in eachinterface(dg, cache)
# Get neighboring element ids
Expand Down Expand Up @@ -93,6 +93,25 @@ end
end

# Calc bounds at physical boundaries
(; boundary_conditions) = semi
calc_bounds_twosided_boundary!(var_min, var_max, variable, u, t,
boundary_conditions,
mesh, equations, dg, cache)

return nothing
end

@inline function calc_bounds_twosided_boundary!(var_min, var_max, variable, u, t,
boundary_conditions::BoundaryConditionPeriodic,
mesh::TreeMesh{2}, equations,
dg, cache)
return nothing
end

@inline function calc_bounds_twosided_boundary!(var_min, var_max, variable, u, t,
boundary_conditions,
mesh::TreeMesh{2}, equations,
dg, cache)
for boundary in eachboundary(dg, cache)
element = cache.boundaries.neighbor_ids[boundary]
orientation = cache.boundaries.orientations[boundary]
Expand Down Expand Up @@ -179,7 +198,7 @@ end
@inline function calc_bounds_onesided_interface!(var_minmax, min_or_max, variable, u, t,
semi, mesh::TreeMesh2D)
_, equations, dg, cache = mesh_equations_solver_cache(semi)
(; boundary_conditions) = semi

# Calc bounds at interfaces and periodic boundaries
for interface in eachinterface(dg, cache)
# Get neighboring element ids
Expand Down Expand Up @@ -214,6 +233,25 @@ end
end

# Calc bounds at physical boundaries
(; boundary_conditions) = semi
calc_bounds_onesided_boundary!(var_minmax, min_or_max, variable, u, t,
boundary_conditions,
mesh, equations, dg, cache)

return nothing
end

@inline function calc_bounds_onesided_boundary!(var_minmax, min_or_max, variable, u, t,
boundary_conditions::BoundaryConditionPeriodic,
mesh::TreeMesh{2}, equations,
dg, cache)
return nothing
end

@inline function calc_bounds_onesided_boundary!(var_minmax, min_or_max, variable, u, t,
boundary_conditions,
mesh::TreeMesh{2}, equations,
dg, cache)
for boundary in eachboundary(dg, cache)
element = cache.boundaries.neighbor_ids[boundary]
orientation = cache.boundaries.orientations[boundary]
Expand Down
131 changes: 131 additions & 0 deletions src/solvers/dgsem_tree/subcell_limiters_3d.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -109,6 +109,72 @@ end
end
end

# Calc bounds at physical boundaries
(; boundary_conditions) = semi
calc_bounds_twosided_boundary!(var_min, var_max, variable, u, t,
boundary_conditions,
mesh, equations, dg, cache)

return nothing
end

@inline function calc_bounds_twosided_boundary!(var_min, var_max, variable,
u, t,
boundary_conditions::BoundaryConditionPeriodic,
mesh::TreeMesh{3}, equations,
dg, cache)
return nothing
end

@inline function calc_bounds_twosided_boundary!(var_min, var_max, variable,
u, t, boundary_conditions,
mesh::TreeMesh{3}, equations,
dg, cache)
for boundary in eachboundary(dg, cache)
element = cache.boundaries.neighbor_ids[boundary]
orientation = cache.boundaries.orientations[boundary]
neighbor_side = cache.boundaries.neighbor_sides[boundary]

for j in eachnode(dg), i in eachnode(dg)
# Define node indices and boundary index based on the orientation and neighbor_side
if neighbor_side == 2 # Element is on the right, boundary on the left
if orientation == 1 # boundary in x-direction
node_index = (1, i, j)
boundary_index = 1
elseif orientation == 2 # boundary in y-direction
node_index = (i, 1, j)
boundary_index = 3
else # orientation == 3 # boundary in z-direction
node_index = (i, j, 1)
boundary_index = 5
end
else # Element is on the left, boundary on the right
if orientation == 1 # boundary in x-direction
node_index = (nnodes(dg), i, j)
boundary_index = 2
elseif orientation == 2 # boundary in y-direction
node_index = (i, nnodes(dg), j)
boundary_index = 4
else # orientation == 3 # boundary in z-direction
node_index = (i, j, nnodes(dg))
boundary_index = 6
end
end
u_inner = get_node_vars(u, equations, dg, node_index..., element)
u_outer = get_boundary_outer_state(u_inner, t,
boundary_conditions[boundary_index],
orientation, boundary_index,
mesh, equations, dg, cache,
node_index..., element)
var_outer = u_outer[variable]

var_min[node_index..., element] = min(var_min[node_index..., element],
var_outer)
var_max[node_index..., element] = max(var_max[node_index..., element],
var_outer)
end
end

return nothing
end

Expand Down Expand Up @@ -211,6 +277,71 @@ end
end
end

# Calc bounds at physical boundaries
(; boundary_conditions) = semi
calc_bounds_onesided_boundary!(var_minmax, min_or_max, variable, u, t,
boundary_conditions,
mesh, equations, dg, cache)

return nothing
end

@inline function calc_bounds_onesided_boundary!(var_minmax, min_or_max, variable,
u, t,
boundary_conditions::BoundaryConditionPeriodic,
mesh::TreeMesh{3}, equations,
dg, cache)
return nothing
end

@inline function calc_bounds_onesided_boundary!(var_minmax, min_or_max, variable,
u, t, boundary_conditions,
mesh::TreeMesh{3}, equations,
dg, cache)
for boundary in eachboundary(dg, cache)
element = cache.boundaries.neighbor_ids[boundary]
orientation = cache.boundaries.orientations[boundary]
neighbor_side = cache.boundaries.neighbor_sides[boundary]

for j in eachnode(dg), i in eachnode(dg)
# Define node indices and boundary index based on the orientation and neighbor_side
if neighbor_side == 2 # Element is on the right, boundary on the left
if orientation == 1 # boundary in x-direction
node_index = (1, i, j)
boundary_index = 1
elseif orientation == 2 # boundary in y-direction
node_index = (i, 1, j)
boundary_index = 3
else # orientation == 3 # boundary in z-direction
node_index = (i, j, 1)
boundary_index = 5
end
else # Element is on the left, boundary on the right
if orientation == 1 # boundary in x-direction
node_index = (nnodes(dg), i, j)
boundary_index = 2
elseif orientation == 2 # boundary in y-direction
node_index = (i, nnodes(dg), j)
boundary_index = 4
else # orientation == 3 # boundary in z-direction
node_index = (i, j, nnodes(dg))
boundary_index = 6
end
end
u_inner = get_node_vars(u, equations, dg, node_index..., element)
u_outer = get_boundary_outer_state(u_inner, t,
boundary_conditions[boundary_index],
orientation, boundary_index,
mesh, equations, dg, cache,
node_index..., element)
var_outer = variable(u_outer, equations)

var_minmax[node_index..., element] = min_or_max(var_minmax[node_index...,
element],
var_outer)
end
end

return nothing
end

Expand Down
Loading