diff --git a/examples/tree_3d_dgsem/elixir_euler_sedov_blast_wave_sc_subcell.jl b/examples/tree_3d_dgsem/elixir_euler_sedov_blast_wave_sc_subcell.jl index ead0b003e81..61452ce3174 100644 --- a/examples/tree_3d_dgsem/elixir_euler_sedov_blast_wave_sc_subcell.jl +++ b/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. diff --git a/src/solvers/dgsem_p4est/subcell_limiters_2d.jl b/src/solvers/dgsem_p4est/subcell_limiters_2d.jl index d42b63789b1..51a9a823d2d 100644 --- a/src/solvers/dgsem_p4est/subcell_limiters_2d.jl +++ b/src/solvers/dgsem_p4est/subcell_limiters_2d.jl @@ -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) @@ -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) @@ -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) @@ -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) diff --git a/src/solvers/dgsem_p4est/subcell_limiters_3d.jl b/src/solvers/dgsem_p4est/subcell_limiters_3d.jl index 4b5d3f1a0ce..fefc3fbbbdc 100644 --- a/src/solvers/dgsem_p4est/subcell_limiters_3d.jl +++ b/src/solvers/dgsem_p4est/subcell_limiters_3d.jl @@ -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) @@ -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) @@ -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) @@ -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) diff --git a/src/solvers/dgsem_structured/subcell_limiters_2d.jl b/src/solvers/dgsem_structured/subcell_limiters_2d.jl index 4be1993bdbe..0623250009e 100644 --- a/src/solvers/dgsem_structured/subcell_limiters_2d.jl +++ b/src/solvers/dgsem_structured/subcell_limiters_2d.jl @@ -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) @@ -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 @@ -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) @@ -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 diff --git a/src/solvers/dgsem_tree/subcell_limiters_2d.jl b/src/solvers/dgsem_tree/subcell_limiters_2d.jl index 2d041ce65ad..e1934e10dfb 100644 --- a/src/solvers/dgsem_tree/subcell_limiters_2d.jl +++ b/src/solvers/dgsem_tree/subcell_limiters_2d.jl @@ -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 @@ -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] @@ -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 @@ -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] diff --git a/src/solvers/dgsem_tree/subcell_limiters_3d.jl b/src/solvers/dgsem_tree/subcell_limiters_3d.jl index 5c8d6c06277..ffe5996d6cd 100644 --- a/src/solvers/dgsem_tree/subcell_limiters_3d.jl +++ b/src/solvers/dgsem_tree/subcell_limiters_3d.jl @@ -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 @@ -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