Add local limiting support for nonperiodic TreeMeshes

examples/tree_3d_dgsem/elixir_euler_sedov_blast_wave_sc_subcell.jl

@@ -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.

src/solvers/dgsem_tree/subcell_limiters_3d.jl

@@ -68,6 +68,7 @@ end
                                                  u, t, semi, mesh::TreeMesh3D,
                                                  equations)
     _, _, dg, cache = mesh_equations_solver_cache(semi)
+    (; boundary_conditions) = semi
 
     # Calc bounds at interfaces and periodic boundaries
     for interface in eachinterface(dg, cache)
@@ -109,6 +110,52 @@ end
         end
     end
 
+    # Calc bounds at physical boundaries
+    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
 
@@ -171,6 +218,7 @@ end
 @inline function calc_bounds_onesided_interface!(var_minmax, min_or_max, variable, u, t,
                                                  semi, mesh::TreeMesh{3})
     _, equations, dg, cache = mesh_equations_solver_cache(semi)
+    (; boundary_conditions) = semi
 
     # Calc bounds at interfaces and periodic boundaries
     for interface in eachinterface(dg, cache)
@@ -211,6 +259,51 @@ end
         end
     end
 
+    # Calc bounds at physical boundaries
+    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