diff --git a/NEWS.md b/NEWS.md index f02da93dbf4..9f74d402707 100644 --- a/NEWS.md +++ b/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 diff --git a/src/solvers/dgsem_p4est/dg.jl b/src/solvers/dgsem_p4est/dg.jl index dddef6999c4..12e13f0d239 100644 --- a/src/solvers/dgsem_p4est/dg.jl +++ b/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") diff --git a/src/solvers/dgsem_p4est/subcell_limiters.jl b/src/solvers/dgsem_p4est/subcell_limiters.jl new file mode 100644 index 00000000000..14b97d91050 --- /dev/null +++ b/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 diff --git a/src/solvers/dgsem_p4est/subcell_limiters_3d.jl b/src/solvers/dgsem_p4est/subcell_limiters_3d.jl index feb4a147136..53490c83f7c 100644 --- a/src/solvers/dgsem_p4est/subcell_limiters_3d.jl +++ b/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 diff --git a/src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl b/src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl index de87fadae85..6ae2765417f 100644 --- a/src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl +++ b/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 diff --git a/src/solvers/dgsem_tree/subcell_limiters.jl b/src/solvers/dgsem_tree/subcell_limiters.jl index 9517ee04859..7a2ecf9df22 100644 --- a/src/solvers/dgsem_tree/subcell_limiters.jl +++ b/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 diff --git a/test/test_p4est_3d.jl b/test/test_p4est_3d.jl index b5908bc8b03..12b2fa92e91 100644 --- a/test/test_p4est_3d.jl +++ b/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=[