Add 3d support for nonperiodic simulations with subcell limiting

NEWS.md

@@ -27,9 +27,10 @@ This enables in particular adaptive mesh refinement for that solver-mesh combina
   Instead, it is possible to reconstruct in custom variables, if functions `cons2recon` and `recon2cons` are provided to
   `VolumeIntegralPureLGLFiniteVolumeO2` and `VolumeIntegralShockCapturingRRG`([#2817]).
 - Add Legendre-Gauss basis for DGSEM and implement solver (`VolumeIntegralWeakForm` and `SurfaceIntegralWeakForm` only) support for conforming 1D & 2D `TreeMesh`es ([#1965]).
-- Extended 3D support for subcell limiting was added ([#2763], [#2846]).
+- Extended 3D support for subcell limiting was added ([#2763], [#2846], [#2844]).
   In the new version, local (minimum and maximum) limiting for conservative variables (using the
   keyword `local_twosided_variables_cons` in `SubcellLimiterIDP()`) with `P4estMesh` is supported.
+  The support was extended to nonperiodic `P4estMesh{3}`es.
   Moreover, initial support for `TreeMesh{3}` including positivity limiting on periodic meshes was added.
 
 ## Changes when updating to v0.15 from v0.14.x

src/solvers/dgsem_p4est/dg.jl

@@ -93,6 +93,8 @@ include("dg_3d.jl")
 include("dg_3d_parabolic.jl")
 include("dg_parallel.jl")
 
+# Subcell limiters
+include("subcell_limiters.jl")
 include("subcell_limiters_2d.jl")
 include("subcell_limiters_3d.jl")
 include("dg_3d_subcell_limiters.jl")

src/solvers/dgsem_p4est/subcell_limiters.jl

@@ -0,0 +1,37 @@
+# By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
+# Since these FMAs can increase the performance of many numerical algorithms,
+# we need to opt-in explicitly.
+# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
+@muladd begin
+#! format: noindent
+
+###############################################################################
+# Auxiliary routine `get_boundary_outer_state` for non-periodic domains
+
+"""
+    get_boundary_outer_state(u_inner, t,
+                             boundary_condition::BoundaryConditionDirichlet,
+                             normal_direction
+                             mesh, equations, dg, cache, indices...)
+For subcell limiting, the calculation of local bounds for non-periodic domains requires the boundary
+outer state. This function returns the boundary value  for [`BoundaryConditionDirichlet`](@ref) at
+time `t` and for node with spatial indices `indices` at the boundary with `normal_direction`.
+
+Should be used together with [`P4estMesh`](@ref).
+
+!!! warning "Experimental implementation"
+    This is an experimental feature and may change in future releases.
+"""
+@inline function get_boundary_outer_state(u_inner, t,
+                                          boundary_condition::BoundaryConditionDirichlet,
+                                          normal_direction,
+                                          mesh::P4estMesh,
+                                          equations, dg, cache, indices...)
+    (; node_coordinates) = cache.elements
+
+    x = get_node_coords(node_coordinates, equations, dg, indices...)
+    u_outer = boundary_condition.boundary_value_function(x, t, equations)
+
+    return u_outer
+end
+end # @muladd

src/solvers/dgsem_p4est/subcell_limiters_3d.jl

@@ -169,6 +169,73 @@ end
     return nothing
 end
 
+@inline function calc_bounds_twosided_boundary!(var_min, var_max, variable, u, t,
+                                                boundary_conditions,
+                                                mesh::P4estMesh{3},
+                                                equations, dg, cache)
+    (; boundary_condition_types, boundary_indices) = boundary_conditions
+    (; contravariant_vectors) = cache.elements
+
+    (; boundaries) = cache
+    index_range = eachnode(dg)
+
+    foreach_enumerate(boundary_condition_types) do (i, boundary_condition)
+        for boundary in boundary_indices[i]
+            element = boundaries.neighbor_ids[boundary]
+            node_indices = boundaries.node_indices[boundary]
+            direction = indices2direction(node_indices)
+
+            i_node_start, i_node_step_i, i_node_step_j = index_to_start_step_3d(node_indices[1],
+                                                                                index_range)
+            j_node_start, j_node_step_i, j_node_step_j = index_to_start_step_3d(node_indices[2],
+                                                                                index_range)
+            k_node_start, k_node_step_i, k_node_step_j = index_to_start_step_3d(node_indices[3],
+                                                                                index_range)
+
+            i_node = i_node_start
+            j_node = j_node_start
+            k_node = k_node_start
+            for j in eachnode(dg)
+                for i in eachnode(dg)
+                    normal_direction = get_normal_direction(direction,
+                                                            contravariant_vectors,
+                                                            i_node, j_node, k_node,
+                                                            element)
+
+                    u_inner = get_node_vars(u, equations, dg, i_node, j_node, k_node,
+                                            element)
+
+                    u_outer = get_boundary_outer_state(u_inner, t, boundary_condition,
+                                                       normal_direction,
+                                                       mesh, equations, dg, cache,
+                                                       i_node, j_node, k_node, element)
+                    var_outer = u_outer[variable]
+
+                    var_min[i_node, j_node, k_node, element] = min(var_min[i_node,
+                                                                           j_node,
+                                                                           k_node,
+                                                                           element],
+                                                                   var_outer)
+                    var_max[i_node, j_node, k_node, element] = max(var_max[i_node,
+                                                                           j_node,
+                                                                           k_node,
+                                                                           element],
+                                                                   var_outer)
+
+                    i_node += i_node_step_i
+                    j_node += j_node_step_i
+                    k_node += k_node_step_i
+                end
+                i_node += i_node_step_j
+                j_node += j_node_step_j
+                k_node += k_node_step_j
+            end
+        end
+    end
+
+    return nothing
+end
+
 @inline function calc_bounds_onesided!(var_minmax, min_or_max, variable,
                                        u::AbstractArray{<:Any, 5}, t, semi)
     mesh, equations, dg, cache = mesh_equations_solver_cache(semi)
@@ -322,6 +389,68 @@ end
     return nothing
 end
 
+@inline function calc_bounds_onesided_boundary!(var_minmax, minmax, variable, u, t,
+                                                boundary_conditions,
+                                                mesh::P4estMesh{3},
+                                                equations, dg, cache)
+    (; boundary_condition_types, boundary_indices) = boundary_conditions
+    (; contravariant_vectors) = cache.elements
+
+    (; boundaries) = cache
+    index_range = eachnode(dg)
+
+    foreach_enumerate(boundary_condition_types) do (i, boundary_condition)
+        for boundary in boundary_indices[i]
+            element = boundaries.neighbor_ids[boundary]
+            node_indices = boundaries.node_indices[boundary]
+            direction = indices2direction(node_indices)
+
+            i_node_start, i_node_step_i, i_node_step_j = index_to_start_step_3d(node_indices[1],
+                                                                                index_range)
+            j_node_start, j_node_step_i, j_node_step_j = index_to_start_step_3d(node_indices[2],
+                                                                                index_range)
+            k_node_start, k_node_step_i, k_node_step_j = index_to_start_step_3d(node_indices[3],
+                                                                                index_range)
+
+            i_node = i_node_start
+            j_node = j_node_start
+            k_node = k_node_start
+            for j in eachnode(dg)
+                for i in eachnode(dg)
+                    normal_direction = get_normal_direction(direction,
+                                                            contravariant_vectors,
+                                                            i_node, j_node, k_node,
+                                                            element)
+
+                    u_inner = get_node_vars(u, equations, dg, i_node, j_node, k_node,
+                                            element)
+
+                    u_outer = get_boundary_outer_state(u_inner, t, boundary_condition,
+                                                       normal_direction,
+                                                       mesh, equations, dg, cache,
+                                                       i_node, j_node, k_node, element)
+                    var_outer = variable(u_outer, equations)
+
+                    var_minmax[i_node, j_node, k_node, element] = minmax(var_minmax[i_node,
+                                                                                    j_node,
+                                                                                    k_node,
+                                                                                    element],
+                                                                         var_outer)
+
+                    i_node += i_node_step_i
+                    j_node += j_node_step_i
+                    k_node += k_node_step_i
+                end
+                i_node += i_node_step_j
+                j_node += j_node_step_j
+                k_node += k_node_step_j
+            end
+        end
+    end
+
+    return nothing
+end
+
 ###############################################################################
 # Local two-sided limiting of conservative variables
 

src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl

@@ -675,32 +675,4 @@ end
 
     return nothing
 end
-
-"""
-    get_boundary_outer_state(u_inner, t,
-                             boundary_condition::BoundaryConditionDirichlet,
-                             orientation_or_normal, direction,
-                             mesh, equations, dg, cache, indices...)
-For subcell limiting, the calculation of local bounds for non-periodic domains requires the boundary
-outer state. This function returns the boundary value  for [`BoundaryConditionDirichlet`](@ref) at
-time `t` and for node with spatial indices `indices` at the boundary with `orientation_or_normal`
-and `direction`.
-
-Should be used together with [`TreeMesh`](@ref) or [`StructuredMesh`](@ref).
-
-!!! warning "Experimental implementation"
-    This is an experimental feature and may change in future releases.
-"""
-@inline function get_boundary_outer_state(u_inner, t,
-                                          boundary_condition::BoundaryConditionDirichlet,
-                                          orientation_or_normal, direction,
-                                          mesh::Union{TreeMesh, StructuredMesh},
-                                          equations, dg, cache, indices...)
-    (; node_coordinates) = cache.elements
-
-    x = get_node_coords(node_coordinates, equations, dg, indices...)
-    u_outer = boundary_condition.boundary_value_function(x, t, equations)
-
-    return u_outer
-end
 end # @muladd

src/solvers/dgsem_tree/subcell_limiters.jl

@@ -447,4 +447,35 @@ end
 @inline function final_check_nonnegative_newton_idp(bound, goal, newton_abstol)
     return (goal <= eps()) && (goal > -max(newton_abstol, abs(bound) * newton_abstol))
 end
+
+###############################################################################
+# Auxiliary routine `get_boundary_outer_state` for non-periodic domains
+
+"""
+    get_boundary_outer_state(u_inner, t,
+                             boundary_condition::BoundaryConditionDirichlet,
+                             orientation_or_normal, direction,
+                             mesh, equations, dg, cache, indices...)
+For subcell limiting, the calculation of local bounds for non-periodic domains requires the boundary
+outer state. This function returns the boundary value  for [`BoundaryConditionDirichlet`](@ref) at
+time `t` and for node with spatial indices `indices` at the boundary with `orientation_or_normal`
+and `direction`.
+
+Should be used together with [`TreeMesh`](@ref) or [`StructuredMesh`](@ref).
+
+!!! warning "Experimental implementation"
+    This is an experimental feature and may change in future releases.
+"""
+@inline function get_boundary_outer_state(u_inner, t,
+                                          boundary_condition::BoundaryConditionDirichlet,
+                                          orientation_or_normal, direction,
+                                          mesh::Union{TreeMesh, StructuredMesh},
+                                          equations, dg, cache, indices...)
+    (; node_coordinates) = cache.elements
+
+    x = get_node_coords(node_coordinates, equations, dg, indices...)
+    u_outer = boundary_condition.boundary_value_function(x, t, equations)
+
+    return u_outer
+end
 end # @muladd

test/test_p4est_3d.jl

@@ -427,6 +427,38 @@ end
     @test_allocations(Trixi.rhs!, semi, sol, 15_000)
 end
 
+@trixi_testset "elixir_euler_sedov_sc_subcell.jl (local bounds, nonperiodic)" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_sedov_sc_subcell.jl"),
+                        local_twosided_variables_cons=["rho"],
+                        local_onesided_variables_nonlinear=[(entropy_guermond_etal,
+                                                             min)],
+                        max_iterations_newton=50,
+                        periodicity=false,
+                        boundary_conditions=BoundaryConditionDirichlet(initial_condition),
+                        l2=[
+                            0.16504564013491585,
+                            0.06461384162458203,
+                            0.06461384162461223,
+                            0.06461384162461678,
+                            0.36193245790622036
+                        ],
+                        linf=[
+                            0.9138327077620716,
+                            0.5707102472596818,
+                            0.5707102472739252,
+                            0.5707102472781822,
+                            4.777595503303726
+                        ],
+                        tspan=(0.0, 0.3))
+    # Ensure that we do not have excessive memory allocations
+    # (e.g., from type instabilities)
+    # Larger values for allowed allocations due to usage of custom
+    # integrator which are not *recorded* for the methods from
+    # OrdinaryDiffEq.jl
+    # Corresponding issue: https://github.com/trixi-framework/Trixi.jl/issues/1877
+    @test_allocations(Trixi.rhs!, semi, sol, 15_000)
+end
+
 @trixi_testset "elixir_euler_sedov.jl (HLLE)" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_sedov.jl"),
                         l2=[