From cde9904f769a5894d4a1c2c9a5d117aa33b7f32d Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Wed, 4 Mar 2026 10:01:17 +0100 Subject: [PATCH 01/46] Start --- src/solvers/dgsem_p4est/containers_2d.jl | 4 +- src/solvers/dgsem_p4est/dg_2d.jl | 158 ++++++++++++++++++ src/solvers/dgsem_structured/containers_2d.jl | 2 +- 3 files changed, 161 insertions(+), 3 deletions(-) diff --git a/src/solvers/dgsem_p4est/containers_2d.jl b/src/solvers/dgsem_p4est/containers_2d.jl index 0cd04609d4a..9f6d31f7ed0 100644 --- a/src/solvers/dgsem_p4est/containers_2d.jl +++ b/src/solvers/dgsem_p4est/containers_2d.jl @@ -8,7 +8,7 @@ # Initialize data structures in element container function init_elements!(elements, mesh::Union{P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}}, - basis::LobattoLegendreBasis) + basis::AbstractBasisSBP) @unpack node_coordinates, jacobian_matrix, contravariant_vectors, inverse_jacobian = elements @@ -29,7 +29,7 @@ end function calc_node_coordinates!(node_coordinates, mesh::Union{P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}}, - basis::LobattoLegendreBasis) + basis::AbstractBasisSBP) # Hanging nodes will cause holes in the mesh if its polydeg is higher # than the polydeg of the solver. @assert length(basis.nodes)>=length(mesh.nodes) "The solver can't have a lower polydeg than the mesh" diff --git a/src/solvers/dgsem_p4est/dg_2d.jl b/src/solvers/dgsem_p4est/dg_2d.jl index 846d5093c6d..4ebab7e130c 100644 --- a/src/solvers/dgsem_p4est/dg_2d.jl +++ b/src/solvers/dgsem_p4est/dg_2d.jl @@ -118,6 +118,107 @@ function prolong2interfaces!(cache, u, return nothing end +function prolong2interfaces!(cache, u, + mesh::Union{P4estMesh{2}, P4estMeshView{2}}, + equations, dg::DGSEM{<:GaussLegendreBasis}) + @unpack interfaces = cache + @unpack boundary_interpolation = dg.basis + index_range = eachnode(dg) + + @threaded for interface in eachinterface(dg, cache) + # Interpolate solution data from the primary element to the interface. + primary_element = interfaces.neighbor_ids[1, interface] + primary_indices = interfaces.node_indices[1, interface] + + i_primary_start, i_primary_step = index_to_start_step_2d(primary_indices[1], + index_range) + j_primary_start, j_primary_step = index_to_start_step_2d(primary_indices[2], + index_range) + # The index direction is identified based on `{i,j}_{primary, secondary}_step`. + # For step = 0, the direction indentified by this index is normal to the face. + # For step != 0 (1 or -1), the direction identified by this index is tangential to the face. + + # Note that in the current implementation, the interface will be + # "aligned at the primary element", i.e., the index of the primary side + # will always run forwards. + + i_primary = i_primary_start + j_primary = j_primary_start + for i in eachnode(dg) + for v in eachvariable(equations) + u_primary = zero(eltype(interfaces.u)) + if i_primary_step == 0 + # i is the normal direction (constant), j varies along the surface + # => Interpolate in first/normal direction + interp_side = (primary_indices[1] === :begin) ? 1 : 2 + # Figure out if we go forward or backward (i.e., query element orientation) + # to know if we need to use `boundary_interpolation[:, 1]` or `boundary_interpolation[:, 2]`. + # If the first index of the element is `:begin`, this face corresponds to left reference element side (-1). + # If the first index of the element is `:end`, this face corresponds to right reference element side (+1). + for ii in eachnode(dg) + u_primary = (u_primary + + u[v, ii, j_primary, primary_element] * + boundary_interpolation[ii, interp_side]) + end + else # j_primary_step == 0 + # j is the normal direction (constant), i varies along the surface + # => Interpolate in second/normal direction + interp_side = (primary_indices[2] === :begin) ? 1 : 2 + for jj in eachnode(dg) + u_primary = (u_primary + + u[v, i_primary, jj, primary_element] * + boundary_interpolation[jj, interp_side]) + end + end + interfaces.u[1, v, i, interface] = u_primary + end + i_primary += i_primary_step + j_primary += j_primary_step + end + + # Interpolate solution data from the secondary element to the interface. + secondary_element = interfaces.neighbor_ids[2, interface] + secondary_indices = interfaces.node_indices[2, interface] + + i_secondary_start, i_secondary_step = index_to_start_step_2d(secondary_indices[1], + index_range) + j_secondary_start, j_secondary_step = index_to_start_step_2d(secondary_indices[2], + index_range) + + i_secondary = i_secondary_start + j_secondary = j_secondary_start + for i in eachnode(dg) + for v in eachvariable(equations) + u_secondary = zero(eltype(interfaces.u)) + if i_secondary_step == 0 + # i is the normal direction (constant), j varies along the surface + # => Interpolate in first/normal direction + interp_side = (secondary_indices[1] === :begin) ? 1 : 2 + for ii in eachnode(dg) + u_secondary = (u_secondary + + u[v, ii, j_secondary, secondary_element] * + boundary_interpolation[ii, interp_side]) + end + else # j_secondary_step == 0 + # j is the normal direction (constant), i varies along the surface + # => Interpolate in second/normal direction + interp_side = (secondary_indices[2] === :begin) ? 1 : 2 + for jj in eachnode(dg) + u_secondary = (u_secondary + + u[v, i_secondary, jj, secondary_element] * + boundary_interpolation[jj, interp_side]) + end + end + interfaces.u[2, v, i, interface] = u_secondary + end + i_secondary += i_secondary_step + j_secondary += j_secondary_step + end + end + + return nothing +end + function calc_interface_flux!(surface_flux_values, mesh::Union{P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}}, @@ -851,6 +952,63 @@ function calc_surface_integral!(du, u, return nothing end +# Specialized variant for `GaussLegendreBasis` that distributes the surface flux +# contribution to all volume nodes. +# Note that all interface fluxes are for `P4estMesh` stored with outward-pointing normals: +# The secondary element's contribution is negated in `calc_interface_flux!` above, +# so all four surface directions are added with a positive sign, matching the existing LGL version above. +# The missing "-" sign is taken care of by `apply_jacobian!`. +function calc_surface_integral!(du, u, + mesh::Union{P4estMesh{2}, P4estMeshView{2}, + T8codeMesh{2}}, + equations, surface_integral::SurfaceIntegralWeakForm, + dg::DGSEM{<:GaussLegendreBasis}, cache) + @unpack boundary_interpolation_inverse_weights = dg.basis + @unpack surface_flux_values = cache.elements + + # We also use explicit assignments instead of `+=` to let `@muladd` turn these + # into FMAs (see comment at the top of the file). + @threaded for element in eachelement(dg, cache) + for l in eachnode(dg) + for v in eachvariable(equations) + # Aliases for the surface flux values on the -x and +x faces + surface_flux_minus_x = surface_flux_values[v, l, 1, element] + surface_flux_plus_x = surface_flux_values[v, l, 2, element] + for ii in eachnode(dg) + # surface at -x: distribute using left boundary interpolation weights + du[v, ii, l, element] = (du[v, ii, l, element] + + surface_flux_minus_x * + boundary_interpolation_inverse_weights[ii, + 1]) + # surface at +x: distribute using right boundary interpolation weights + du[v, ii, l, element] = (du[v, ii, l, element] + + surface_flux_plus_x * + boundary_interpolation_inverse_weights[ii, + 2]) + end + + # Aliases for the surface flux values on the -y and +y faces + surface_flux_minus_y = surface_flux_values[v, l, 3, element] + surface_flux_plus_y = surface_flux_values[v, l, 4, element] + for jj in eachnode(dg) + # surface at -y: distribute using left boundary interpolation weights + du[v, l, jj, element] = (du[v, l, jj, element] + + surface_flux_minus_y * + boundary_interpolation_inverse_weights[jj, + 1]) + # surface at +y: distribute using right boundary interpolation weights + du[v, l, jj, element] = (du[v, l, jj, element] + + surface_flux_plus_y * + boundary_interpolation_inverse_weights[jj, + 2]) + end + end + end + end + + return nothing +end + # Call this for coupled P4estMeshView simulations. # The coupling calculations (especially boundary conditions) require data from the parent mesh, which is why # the additional variable u_parent is needed, compared to non-coupled systems. diff --git a/src/solvers/dgsem_structured/containers_2d.jl b/src/solvers/dgsem_structured/containers_2d.jl index 71151a785d2..03241705266 100644 --- a/src/solvers/dgsem_structured/containers_2d.jl +++ b/src/solvers/dgsem_structured/containers_2d.jl @@ -60,7 +60,7 @@ end # Calculate Jacobian matrix of the mapping from the reference element to the element in the physical domain function calc_jacobian_matrix!(jacobian_matrix, element, node_coordinates::AbstractArray{<:Any, 4}, - basis::LobattoLegendreBasis) + basis::AbstractBasisSBP) @unpack derivative_matrix = basis # The code below is equivalent to the following matrix multiplications, which From ad693080c53a771ad5a82d12e6151760797627ef Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Wed, 4 Mar 2026 10:02:19 +0100 Subject: [PATCH 02/46] Revert "Start" This reverts commit cde9904f769a5894d4a1c2c9a5d117aa33b7f32d. --- src/solvers/dgsem_p4est/containers_2d.jl | 4 +- src/solvers/dgsem_p4est/dg_2d.jl | 158 ------------------ src/solvers/dgsem_structured/containers_2d.jl | 2 +- 3 files changed, 3 insertions(+), 161 deletions(-) diff --git a/src/solvers/dgsem_p4est/containers_2d.jl b/src/solvers/dgsem_p4est/containers_2d.jl index 9f6d31f7ed0..0cd04609d4a 100644 --- a/src/solvers/dgsem_p4est/containers_2d.jl +++ b/src/solvers/dgsem_p4est/containers_2d.jl @@ -8,7 +8,7 @@ # Initialize data structures in element container function init_elements!(elements, mesh::Union{P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}}, - basis::AbstractBasisSBP) + basis::LobattoLegendreBasis) @unpack node_coordinates, jacobian_matrix, contravariant_vectors, inverse_jacobian = elements @@ -29,7 +29,7 @@ end function calc_node_coordinates!(node_coordinates, mesh::Union{P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}}, - basis::AbstractBasisSBP) + basis::LobattoLegendreBasis) # Hanging nodes will cause holes in the mesh if its polydeg is higher # than the polydeg of the solver. @assert length(basis.nodes)>=length(mesh.nodes) "The solver can't have a lower polydeg than the mesh" diff --git a/src/solvers/dgsem_p4est/dg_2d.jl b/src/solvers/dgsem_p4est/dg_2d.jl index 4ebab7e130c..846d5093c6d 100644 --- a/src/solvers/dgsem_p4est/dg_2d.jl +++ b/src/solvers/dgsem_p4est/dg_2d.jl @@ -118,107 +118,6 @@ function prolong2interfaces!(cache, u, return nothing end -function prolong2interfaces!(cache, u, - mesh::Union{P4estMesh{2}, P4estMeshView{2}}, - equations, dg::DGSEM{<:GaussLegendreBasis}) - @unpack interfaces = cache - @unpack boundary_interpolation = dg.basis - index_range = eachnode(dg) - - @threaded for interface in eachinterface(dg, cache) - # Interpolate solution data from the primary element to the interface. - primary_element = interfaces.neighbor_ids[1, interface] - primary_indices = interfaces.node_indices[1, interface] - - i_primary_start, i_primary_step = index_to_start_step_2d(primary_indices[1], - index_range) - j_primary_start, j_primary_step = index_to_start_step_2d(primary_indices[2], - index_range) - # The index direction is identified based on `{i,j}_{primary, secondary}_step`. - # For step = 0, the direction indentified by this index is normal to the face. - # For step != 0 (1 or -1), the direction identified by this index is tangential to the face. - - # Note that in the current implementation, the interface will be - # "aligned at the primary element", i.e., the index of the primary side - # will always run forwards. - - i_primary = i_primary_start - j_primary = j_primary_start - for i in eachnode(dg) - for v in eachvariable(equations) - u_primary = zero(eltype(interfaces.u)) - if i_primary_step == 0 - # i is the normal direction (constant), j varies along the surface - # => Interpolate in first/normal direction - interp_side = (primary_indices[1] === :begin) ? 1 : 2 - # Figure out if we go forward or backward (i.e., query element orientation) - # to know if we need to use `boundary_interpolation[:, 1]` or `boundary_interpolation[:, 2]`. - # If the first index of the element is `:begin`, this face corresponds to left reference element side (-1). - # If the first index of the element is `:end`, this face corresponds to right reference element side (+1). - for ii in eachnode(dg) - u_primary = (u_primary + - u[v, ii, j_primary, primary_element] * - boundary_interpolation[ii, interp_side]) - end - else # j_primary_step == 0 - # j is the normal direction (constant), i varies along the surface - # => Interpolate in second/normal direction - interp_side = (primary_indices[2] === :begin) ? 1 : 2 - for jj in eachnode(dg) - u_primary = (u_primary + - u[v, i_primary, jj, primary_element] * - boundary_interpolation[jj, interp_side]) - end - end - interfaces.u[1, v, i, interface] = u_primary - end - i_primary += i_primary_step - j_primary += j_primary_step - end - - # Interpolate solution data from the secondary element to the interface. - secondary_element = interfaces.neighbor_ids[2, interface] - secondary_indices = interfaces.node_indices[2, interface] - - i_secondary_start, i_secondary_step = index_to_start_step_2d(secondary_indices[1], - index_range) - j_secondary_start, j_secondary_step = index_to_start_step_2d(secondary_indices[2], - index_range) - - i_secondary = i_secondary_start - j_secondary = j_secondary_start - for i in eachnode(dg) - for v in eachvariable(equations) - u_secondary = zero(eltype(interfaces.u)) - if i_secondary_step == 0 - # i is the normal direction (constant), j varies along the surface - # => Interpolate in first/normal direction - interp_side = (secondary_indices[1] === :begin) ? 1 : 2 - for ii in eachnode(dg) - u_secondary = (u_secondary + - u[v, ii, j_secondary, secondary_element] * - boundary_interpolation[ii, interp_side]) - end - else # j_secondary_step == 0 - # j is the normal direction (constant), i varies along the surface - # => Interpolate in second/normal direction - interp_side = (secondary_indices[2] === :begin) ? 1 : 2 - for jj in eachnode(dg) - u_secondary = (u_secondary + - u[v, i_secondary, jj, secondary_element] * - boundary_interpolation[jj, interp_side]) - end - end - interfaces.u[2, v, i, interface] = u_secondary - end - i_secondary += i_secondary_step - j_secondary += j_secondary_step - end - end - - return nothing -end - function calc_interface_flux!(surface_flux_values, mesh::Union{P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}}, @@ -952,63 +851,6 @@ function calc_surface_integral!(du, u, return nothing end -# Specialized variant for `GaussLegendreBasis` that distributes the surface flux -# contribution to all volume nodes. -# Note that all interface fluxes are for `P4estMesh` stored with outward-pointing normals: -# The secondary element's contribution is negated in `calc_interface_flux!` above, -# so all four surface directions are added with a positive sign, matching the existing LGL version above. -# The missing "-" sign is taken care of by `apply_jacobian!`. -function calc_surface_integral!(du, u, - mesh::Union{P4estMesh{2}, P4estMeshView{2}, - T8codeMesh{2}}, - equations, surface_integral::SurfaceIntegralWeakForm, - dg::DGSEM{<:GaussLegendreBasis}, cache) - @unpack boundary_interpolation_inverse_weights = dg.basis - @unpack surface_flux_values = cache.elements - - # We also use explicit assignments instead of `+=` to let `@muladd` turn these - # into FMAs (see comment at the top of the file). - @threaded for element in eachelement(dg, cache) - for l in eachnode(dg) - for v in eachvariable(equations) - # Aliases for the surface flux values on the -x and +x faces - surface_flux_minus_x = surface_flux_values[v, l, 1, element] - surface_flux_plus_x = surface_flux_values[v, l, 2, element] - for ii in eachnode(dg) - # surface at -x: distribute using left boundary interpolation weights - du[v, ii, l, element] = (du[v, ii, l, element] + - surface_flux_minus_x * - boundary_interpolation_inverse_weights[ii, - 1]) - # surface at +x: distribute using right boundary interpolation weights - du[v, ii, l, element] = (du[v, ii, l, element] + - surface_flux_plus_x * - boundary_interpolation_inverse_weights[ii, - 2]) - end - - # Aliases for the surface flux values on the -y and +y faces - surface_flux_minus_y = surface_flux_values[v, l, 3, element] - surface_flux_plus_y = surface_flux_values[v, l, 4, element] - for jj in eachnode(dg) - # surface at -y: distribute using left boundary interpolation weights - du[v, l, jj, element] = (du[v, l, jj, element] + - surface_flux_minus_y * - boundary_interpolation_inverse_weights[jj, - 1]) - # surface at +y: distribute using right boundary interpolation weights - du[v, l, jj, element] = (du[v, l, jj, element] + - surface_flux_plus_y * - boundary_interpolation_inverse_weights[jj, - 2]) - end - end - end - end - - return nothing -end - # Call this for coupled P4estMeshView simulations. # The coupling calculations (especially boundary conditions) require data from the parent mesh, which is why # the additional variable u_parent is needed, compared to non-coupled systems. diff --git a/src/solvers/dgsem_structured/containers_2d.jl b/src/solvers/dgsem_structured/containers_2d.jl index 03241705266..71151a785d2 100644 --- a/src/solvers/dgsem_structured/containers_2d.jl +++ b/src/solvers/dgsem_structured/containers_2d.jl @@ -60,7 +60,7 @@ end # Calculate Jacobian matrix of the mapping from the reference element to the element in the physical domain function calc_jacobian_matrix!(jacobian_matrix, element, node_coordinates::AbstractArray{<:Any, 4}, - basis::AbstractBasisSBP) + basis::LobattoLegendreBasis) @unpack derivative_matrix = basis # The code below is equivalent to the following matrix multiplications, which From 8baf1704c8950649c19de1244f930adc16401048 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Wed, 4 Mar 2026 10:10:24 +0100 Subject: [PATCH 03/46] Start --- src/solvers/dgsem_p4est/containers_2d.jl | 4 +- src/solvers/dgsem_p4est/dg_2d.jl | 160 +++++++++++++++++- src/solvers/dgsem_structured/containers_2d.jl | 2 +- test/test_p4est_2d.jl | 12 ++ 4 files changed, 174 insertions(+), 4 deletions(-) diff --git a/src/solvers/dgsem_p4est/containers_2d.jl b/src/solvers/dgsem_p4est/containers_2d.jl index 0cd04609d4a..9f6d31f7ed0 100644 --- a/src/solvers/dgsem_p4est/containers_2d.jl +++ b/src/solvers/dgsem_p4est/containers_2d.jl @@ -8,7 +8,7 @@ # Initialize data structures in element container function init_elements!(elements, mesh::Union{P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}}, - basis::LobattoLegendreBasis) + basis::AbstractBasisSBP) @unpack node_coordinates, jacobian_matrix, contravariant_vectors, inverse_jacobian = elements @@ -29,7 +29,7 @@ end function calc_node_coordinates!(node_coordinates, mesh::Union{P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}}, - basis::LobattoLegendreBasis) + basis::AbstractBasisSBP) # Hanging nodes will cause holes in the mesh if its polydeg is higher # than the polydeg of the solver. @assert length(basis.nodes)>=length(mesh.nodes) "The solver can't have a lower polydeg than the mesh" diff --git a/src/solvers/dgsem_p4est/dg_2d.jl b/src/solvers/dgsem_p4est/dg_2d.jl index 846d5093c6d..c4afd69b6b9 100644 --- a/src/solvers/dgsem_p4est/dg_2d.jl +++ b/src/solvers/dgsem_p4est/dg_2d.jl @@ -118,6 +118,107 @@ function prolong2interfaces!(cache, u, return nothing end +function prolong2interfaces!(cache, u, + mesh::Union{P4estMesh{2}, P4estMeshView{2}}, + equations, dg::DGSEM{<:GaussLegendreBasis}) + @unpack interfaces = cache + @unpack boundary_interpolation = dg.basis + index_range = eachnode(dg) + + @threaded for interface in eachinterface(dg, cache) + # Interpolate solution data from the primary element to the interface. + primary_element = interfaces.neighbor_ids[1, interface] + primary_indices = interfaces.node_indices[1, interface] + + i_primary_start, i_primary_step = index_to_start_step_2d(primary_indices[1], + index_range) + j_primary_start, j_primary_step = index_to_start_step_2d(primary_indices[2], + index_range) + # The index direction is identified based on `{i,j}_{primary, secondary}_step`. + # For step = 0, the direction indentified by this index is normal to the face. + # For step != 0 (1 or -1), the direction identified by this index is tangential to the face. + + # Note that in the current implementation, the interface will be + # "aligned at the primary element", i.e., the index of the primary side + # will always run forwards. + + i_primary = i_primary_start + j_primary = j_primary_start + for i in eachnode(dg) + for v in eachvariable(equations) + u_primary = zero(eltype(interfaces.u)) + if i_primary_step == 0 + # i is the normal direction (constant), j varies along the surface + # => Interpolate in first/normal direction + interp_side = (primary_indices[1] === :begin) ? 1 : 2 + # Figure out if we go forward or backward (i.e., query element orientation) + # to know if we need to use `boundary_interpolation[:, 1]` or `boundary_interpolation[:, 2]`. + # If the first index of the element is `:begin`, this face corresponds to left reference element side (-1). + # If the first index of the element is `:end`, this face corresponds to right reference element side (+1). + for ii in eachnode(dg) + u_primary = (u_primary + + u[v, ii, j_primary, primary_element] * + boundary_interpolation[ii, interp_side]) + end + else # j_primary_step == 0 + # j is the normal direction (constant), i varies along the surface + # => Interpolate in second/normal direction + interp_side = (primary_indices[2] === :begin) ? 1 : 2 + for jj in eachnode(dg) + u_primary = (u_primary + + u[v, i_primary, jj, primary_element] * + boundary_interpolation[jj, interp_side]) + end + end + interfaces.u[1, v, i, interface] = u_primary + end + i_primary += i_primary_step + j_primary += j_primary_step + end + + # Interpolate solution data from the secondary element to the interface. + secondary_element = interfaces.neighbor_ids[2, interface] + secondary_indices = interfaces.node_indices[2, interface] + + i_secondary_start, i_secondary_step = index_to_start_step_2d(secondary_indices[1], + index_range) + j_secondary_start, j_secondary_step = index_to_start_step_2d(secondary_indices[2], + index_range) + + i_secondary = i_secondary_start + j_secondary = j_secondary_start + for i in eachnode(dg) + for v in eachvariable(equations) + u_secondary = zero(eltype(interfaces.u)) + if i_secondary_step == 0 + # i is the normal direction (constant), j varies along the surface + # => Interpolate in first/normal direction + interp_side = (secondary_indices[1] === :begin) ? 1 : 2 + for ii in eachnode(dg) + u_secondary = (u_secondary + + u[v, ii, j_secondary, secondary_element] * + boundary_interpolation[ii, interp_side]) + end + else # j_secondary_step == 0 + # j is the normal direction (constant), i varies along the surface + # => Interpolate in second/normal direction + interp_side = (secondary_indices[2] === :begin) ? 1 : 2 + for jj in eachnode(dg) + u_secondary = (u_secondary + + u[v, i_secondary, jj, secondary_element] * + boundary_interpolation[jj, interp_side]) + end + end + interfaces.u[2, v, i, interface] = u_secondary + end + i_secondary += i_secondary_step + j_secondary += j_secondary_step + end + end + + return nothing +end + function calc_interface_flux!(surface_flux_values, mesh::Union{P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}}, @@ -851,6 +952,63 @@ function calc_surface_integral!(du, u, return nothing end +# Specialized variant for `GaussLegendreBasis` that distributes the surface flux +# contribution to all volume nodes. +# Note that all interface fluxes are for `P4estMesh` stored with outward-pointing normals: +# The secondary element's contribution is negated in `calc_interface_flux!` above, +# so all four surface directions are added with a positive sign, matching the existing LGL version above. +# The missing "-" sign is taken care of by `apply_jacobian!`. +function calc_surface_integral!(du, u, + mesh::Union{P4estMesh{2}, P4estMeshView{2}, + T8codeMesh{2}}, + equations, surface_integral::SurfaceIntegralWeakForm, + dg::DGSEM{<:GaussLegendreBasis}, cache) + @unpack boundary_interpolation_inverse_weights = dg.basis + @unpack surface_flux_values = cache.elements + + # We also use explicit assignments instead of `+=` to let `@muladd` turn these + # into FMAs (see comment at the top of the file). + @threaded for element in eachelement(dg, cache) + for l in eachnode(dg) + for v in eachvariable(equations) + # Aliases for the surface flux values on the -x and +x faces + surface_flux_minus_x = surface_flux_values[v, l, 1, element] + surface_flux_plus_x = surface_flux_values[v, l, 2, element] + for ii in eachnode(dg) + # surface at -x: distribute using left boundary interpolation weights + du[v, ii, l, element] = (du[v, ii, l, element] + + surface_flux_minus_x * + boundary_interpolation_inverse_weights[ii, + 1]) + # surface at +x: distribute using right boundary interpolation weights + du[v, ii, l, element] = (du[v, ii, l, element] + + surface_flux_plus_x * + boundary_interpolation_inverse_weights[ii, + 2]) + end + + # Aliases for the surface flux values on the -y and +y faces + surface_flux_minus_y = surface_flux_values[v, l, 3, element] + surface_flux_plus_y = surface_flux_values[v, l, 4, element] + for jj in eachnode(dg) + # surface at -y: distribute using left boundary interpolation weights + du[v, l, jj, element] = (du[v, l, jj, element] + + surface_flux_minus_y * + boundary_interpolation_inverse_weights[jj, + 1]) + # surface at +y: distribute using right boundary interpolation weights + du[v, l, jj, element] = (du[v, l, jj, element] + + surface_flux_plus_y * + boundary_interpolation_inverse_weights[jj, + 2]) + end + end + end + end + + return nothing +end + # Call this for coupled P4estMeshView simulations. # The coupling calculations (especially boundary conditions) require data from the parent mesh, which is why # the additional variable u_parent is needed, compared to non-coupled systems. @@ -922,4 +1080,4 @@ function rhs!(du, u, t, u_parent, semis, return nothing end -end # @muladd +end # @muladd \ No newline at end of file diff --git a/src/solvers/dgsem_structured/containers_2d.jl b/src/solvers/dgsem_structured/containers_2d.jl index 71151a785d2..03241705266 100644 --- a/src/solvers/dgsem_structured/containers_2d.jl +++ b/src/solvers/dgsem_structured/containers_2d.jl @@ -60,7 +60,7 @@ end # Calculate Jacobian matrix of the mapping from the reference element to the element in the physical domain function calc_jacobian_matrix!(jacobian_matrix, element, node_coordinates::AbstractArray{<:Any, 4}, - basis::LobattoLegendreBasis) + basis::AbstractBasisSBP) @unpack derivative_matrix = basis # The code below is equivalent to the following matrix multiplications, which diff --git a/test/test_p4est_2d.jl b/test/test_p4est_2d.jl index 4ebd2599314..b4f78cc8651 100644 --- a/test/test_p4est_2d.jl +++ b/test/test_p4est_2d.jl @@ -30,6 +30,18 @@ isdir(outdir) && rm(outdir, recursive = true) @test real(semi32.mesh) == Float64 end +@trixi_testset "elixir_advection_basic.jl (Gauss-Legendre)" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_basic.jl"), + solver = DGSEM(polydeg = 3, surface_flux = flux_lax_friedrichs, + basis_type = GaussLegendreBasis), + cfl = 0.8, + l2=[8.311947673061856e-6], + linf=[6.627000273229378e-5]) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + @test_allocations(Trixi.rhs!, semi, sol, 1000) +end + @trixi_testset "elixir_advection_basic.jl (Float32)" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_basic_gpu.jl"), # Expected errors are exactly the same as with TreeMesh! From a20bdb2eca99757aed38dd044eacf1709d58052a Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Wed, 4 Mar 2026 10:16:51 +0100 Subject: [PATCH 04/46] test --- src/solvers/dgsem_p4est/dg_2d.jl | 11 +++++------ test/test_p4est_2d.jl | 9 ++++----- 2 files changed, 9 insertions(+), 11 deletions(-) diff --git a/src/solvers/dgsem_p4est/dg_2d.jl b/src/solvers/dgsem_p4est/dg_2d.jl index c4afd69b6b9..6a272ce6423 100644 --- a/src/solvers/dgsem_p4est/dg_2d.jl +++ b/src/solvers/dgsem_p4est/dg_2d.jl @@ -952,12 +952,6 @@ function calc_surface_integral!(du, u, return nothing end -# Specialized variant for `GaussLegendreBasis` that distributes the surface flux -# contribution to all volume nodes. -# Note that all interface fluxes are for `P4estMesh` stored with outward-pointing normals: -# The secondary element's contribution is negated in `calc_interface_flux!` above, -# so all four surface directions are added with a positive sign, matching the existing LGL version above. -# The missing "-" sign is taken care of by `apply_jacobian!`. function calc_surface_integral!(du, u, mesh::Union{P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}}, @@ -966,6 +960,11 @@ function calc_surface_integral!(du, u, @unpack boundary_interpolation_inverse_weights = dg.basis @unpack surface_flux_values = cache.elements + # Note that all fluxes have been computed with outward-pointing normal vectors. + # This computes the **negative** surface integral contribution, + # i.e., M^{-1} * boundary_interpolation^T + # and the missing "-" is taken care of by `apply_jacobian!`. + # # We also use explicit assignments instead of `+=` to let `@muladd` turn these # into FMAs (see comment at the top of the file). @threaded for element in eachelement(dg, cache) diff --git a/test/test_p4est_2d.jl b/test/test_p4est_2d.jl index b4f78cc8651..b030dfe4fe4 100644 --- a/test/test_p4est_2d.jl +++ b/test/test_p4est_2d.jl @@ -32,11 +32,10 @@ end @trixi_testset "elixir_advection_basic.jl (Gauss-Legendre)" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_basic.jl"), - solver = DGSEM(polydeg = 3, surface_flux = flux_lax_friedrichs, - basis_type = GaussLegendreBasis), - cfl = 0.8, - l2=[8.311947673061856e-6], - linf=[6.627000273229378e-5]) + solver=DGSEM(polydeg = 3, surface_flux = flux_lax_friedrichs, + basis_type = GaussLegendreBasis), + cfl=0.8, + l2=[3.721398353159235e-6], linf=[1.8621131703255855e-5]) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) @test_allocations(Trixi.rhs!, semi, sol, 1000) From 9a8fc4dd10e43d76f77282b4ae4f15d8c2ab9736 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Wed, 4 Mar 2026 10:17:42 +0100 Subject: [PATCH 05/46] fmt --- src/solvers/dgsem_p4est/dg_2d.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/solvers/dgsem_p4est/dg_2d.jl b/src/solvers/dgsem_p4est/dg_2d.jl index 6a272ce6423..2c8b22d7607 100644 --- a/src/solvers/dgsem_p4est/dg_2d.jl +++ b/src/solvers/dgsem_p4est/dg_2d.jl @@ -1079,4 +1079,4 @@ function rhs!(du, u, t, u_parent, semis, return nothing end -end # @muladd \ No newline at end of file +end # @muladd From 4c4054671d7b4aa2ae271106452747c83c8da3e8 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Wed, 4 Mar 2026 10:24:45 +0100 Subject: [PATCH 06/46] restrict --- src/solvers/dgsem_p4est/dg_2d.jl | 14 ++++++-------- src/solvers/dgsem_tree/dg_1d.jl | 2 +- src/solvers/dgsem_tree/dg_2d.jl | 3 +-- 3 files changed, 8 insertions(+), 11 deletions(-) diff --git a/src/solvers/dgsem_p4est/dg_2d.jl b/src/solvers/dgsem_p4est/dg_2d.jl index 2c8b22d7607..0ebda59b05a 100644 --- a/src/solvers/dgsem_p4est/dg_2d.jl +++ b/src/solvers/dgsem_p4est/dg_2d.jl @@ -953,8 +953,7 @@ function calc_surface_integral!(du, u, end function calc_surface_integral!(du, u, - mesh::Union{P4estMesh{2}, P4estMeshView{2}, - T8codeMesh{2}}, + mesh::Union{P4estMesh{2}, P4estMeshView{2}}, equations, surface_integral::SurfaceIntegralWeakForm, dg::DGSEM{<:GaussLegendreBasis}, cache) @unpack boundary_interpolation_inverse_weights = dg.basis @@ -970,32 +969,31 @@ function calc_surface_integral!(du, u, @threaded for element in eachelement(dg, cache) for l in eachnode(dg) for v in eachvariable(equations) - # Aliases for the surface flux values on the -x and +x faces + # Aliases for repeatedly accessed variables surface_flux_minus_x = surface_flux_values[v, l, 1, element] surface_flux_plus_x = surface_flux_values[v, l, 2, element] for ii in eachnode(dg) - # surface at -x: distribute using left boundary interpolation weights + # surface at -x du[v, ii, l, element] = (du[v, ii, l, element] + surface_flux_minus_x * boundary_interpolation_inverse_weights[ii, 1]) - # surface at +x: distribute using right boundary interpolation weights + # surface at +x du[v, ii, l, element] = (du[v, ii, l, element] + surface_flux_plus_x * boundary_interpolation_inverse_weights[ii, 2]) end - # Aliases for the surface flux values on the -y and +y faces surface_flux_minus_y = surface_flux_values[v, l, 3, element] surface_flux_plus_y = surface_flux_values[v, l, 4, element] for jj in eachnode(dg) - # surface at -y: distribute using left boundary interpolation weights + # surface at -y du[v, l, jj, element] = (du[v, l, jj, element] + surface_flux_minus_y * boundary_interpolation_inverse_weights[jj, 1]) - # surface at +y: distribute using right boundary interpolation weights + # surface at +y du[v, l, jj, element] = (du[v, l, jj, element] + surface_flux_plus_y * boundary_interpolation_inverse_weights[jj, diff --git a/src/solvers/dgsem_tree/dg_1d.jl b/src/solvers/dgsem_tree/dg_1d.jl index 58664b2ea98..cd8e6c3fbda 100644 --- a/src/solvers/dgsem_tree/dg_1d.jl +++ b/src/solvers/dgsem_tree/dg_1d.jl @@ -696,7 +696,7 @@ function calc_surface_integral!(du, u, mesh::Union{TreeMesh{1}, StructuredMesh{1 return nothing end -function calc_surface_integral!(du, u, mesh::Union{TreeMesh{1}, StructuredMesh{1}}, +function calc_surface_integral!(du, u, mesh::TreeMesh{1}, equations, surface_integral::SurfaceIntegralWeakForm, dg::DGSEM{<:GaussLegendreBasis}, cache) @unpack boundary_interpolation_inverse_weights = dg.basis diff --git a/src/solvers/dgsem_tree/dg_2d.jl b/src/solvers/dgsem_tree/dg_2d.jl index a944fdca4b6..bcd6249f154 100644 --- a/src/solvers/dgsem_tree/dg_2d.jl +++ b/src/solvers/dgsem_tree/dg_2d.jl @@ -1249,8 +1249,7 @@ function calc_surface_integral!(du, u, end function calc_surface_integral!(du, u, - mesh::Union{TreeMesh{2}, StructuredMesh{2}, - StructuredMeshView{2}}, + mesh::TreeMesh{2}, equations, surface_integral::SurfaceIntegralWeakForm, dg::DGSEM{<:GaussLegendreBasis}, cache) @unpack boundary_interpolation_inverse_weights = dg.basis From 9deee4788ab2bbba918b57814c45839785c8e929 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Wed, 4 Mar 2026 10:27:23 +0100 Subject: [PATCH 07/46] typo --- src/solvers/dgsem_p4est/dg_2d.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/solvers/dgsem_p4est/dg_2d.jl b/src/solvers/dgsem_p4est/dg_2d.jl index 0ebda59b05a..9b5194ab109 100644 --- a/src/solvers/dgsem_p4est/dg_2d.jl +++ b/src/solvers/dgsem_p4est/dg_2d.jl @@ -135,7 +135,7 @@ function prolong2interfaces!(cache, u, j_primary_start, j_primary_step = index_to_start_step_2d(primary_indices[2], index_range) # The index direction is identified based on `{i,j}_{primary, secondary}_step`. - # For step = 0, the direction indentified by this index is normal to the face. + # For step = 0, the direction identified by this index is normal to the face. # For step != 0 (1 or -1), the direction identified by this index is tangential to the face. # Note that in the current implementation, the interface will be From 883d999cea423ed79f451fdfc2d8ae32ba91c02a Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Mon, 9 Mar 2026 16:01:19 +0100 Subject: [PATCH 08/46] pre-comp --- src/solvers/dgsem_p4est/dg_2d.jl | 39 +++++++++++++++++++------------- 1 file changed, 23 insertions(+), 16 deletions(-) diff --git a/src/solvers/dgsem_p4est/dg_2d.jl b/src/solvers/dgsem_p4est/dg_2d.jl index 9b5194ab109..6b18415316c 100644 --- a/src/solvers/dgsem_p4est/dg_2d.jl +++ b/src/solvers/dgsem_p4est/dg_2d.jl @@ -138,6 +138,16 @@ function prolong2interfaces!(cache, u, # For step = 0, the direction identified by this index is normal to the face. # For step != 0 (1 or -1), the direction identified by this index is tangential to the face. + # Precompute the interpolation indices for the primary element, case: + + # 1) i is the normal direction (constant), j varies along the surface + # => Interpolate in first/normal direction + interp_side_1 = (primary_indices[1] === :begin) ? 1 : 2 + + # 2) j is the normal direction (constant), i varies along the surface + # => Interpolate in second/normal direction + interp_side_2 = (primary_indices[2] === :begin) ? 1 : 2 + # Note that in the current implementation, the interface will be # "aligned at the primary element", i.e., the index of the primary side # will always run forwards. @@ -148,9 +158,6 @@ function prolong2interfaces!(cache, u, for v in eachvariable(equations) u_primary = zero(eltype(interfaces.u)) if i_primary_step == 0 - # i is the normal direction (constant), j varies along the surface - # => Interpolate in first/normal direction - interp_side = (primary_indices[1] === :begin) ? 1 : 2 # Figure out if we go forward or backward (i.e., query element orientation) # to know if we need to use `boundary_interpolation[:, 1]` or `boundary_interpolation[:, 2]`. # If the first index of the element is `:begin`, this face corresponds to left reference element side (-1). @@ -158,16 +165,13 @@ function prolong2interfaces!(cache, u, for ii in eachnode(dg) u_primary = (u_primary + u[v, ii, j_primary, primary_element] * - boundary_interpolation[ii, interp_side]) + boundary_interpolation[ii, interp_side_1]) end else # j_primary_step == 0 - # j is the normal direction (constant), i varies along the surface - # => Interpolate in second/normal direction - interp_side = (primary_indices[2] === :begin) ? 1 : 2 for jj in eachnode(dg) u_primary = (u_primary + u[v, i_primary, jj, primary_element] * - boundary_interpolation[jj, interp_side]) + boundary_interpolation[jj, interp_side_2]) end end interfaces.u[1, v, i, interface] = u_primary @@ -185,28 +189,31 @@ function prolong2interfaces!(cache, u, j_secondary_start, j_secondary_step = index_to_start_step_2d(secondary_indices[2], index_range) + # Precompute the interpolation indices for the secondary element, case: + # 1) i is the normal direction (constant), j varies along the surface + # => Interpolate in first/normal direction + interp_side_1 = (secondary_indices[1] === :begin) ? 1 : 2 + + # 2) j is the normal direction (constant), i varies along the surface + # => Interpolate in second/normal direction + interp_side_2 = (secondary_indices[2] === :begin) ? 1 : 2 + i_secondary = i_secondary_start j_secondary = j_secondary_start for i in eachnode(dg) for v in eachvariable(equations) u_secondary = zero(eltype(interfaces.u)) if i_secondary_step == 0 - # i is the normal direction (constant), j varies along the surface - # => Interpolate in first/normal direction - interp_side = (secondary_indices[1] === :begin) ? 1 : 2 for ii in eachnode(dg) u_secondary = (u_secondary + u[v, ii, j_secondary, secondary_element] * - boundary_interpolation[ii, interp_side]) + boundary_interpolation[ii, interp_side_1]) end else # j_secondary_step == 0 - # j is the normal direction (constant), i varies along the surface - # => Interpolate in second/normal direction - interp_side = (secondary_indices[2] === :begin) ? 1 : 2 for jj in eachnode(dg) u_secondary = (u_secondary + u[v, i_secondary, jj, secondary_element] * - boundary_interpolation[jj, interp_side]) + boundary_interpolation[jj, interp_side_2]) end end interfaces.u[2, v, i, interface] = u_secondary From 89efbbe635a523c7f28fdbac908f2cc2d76480a1 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Mon, 9 Mar 2026 16:08:46 +0100 Subject: [PATCH 09/46] optimize --- src/solvers/dgsem_p4est/dg_2d.jl | 89 +++++++++++++++++--------------- 1 file changed, 46 insertions(+), 43 deletions(-) diff --git a/src/solvers/dgsem_p4est/dg_2d.jl b/src/solvers/dgsem_p4est/dg_2d.jl index 6b18415316c..61a865c4f37 100644 --- a/src/solvers/dgsem_p4est/dg_2d.jl +++ b/src/solvers/dgsem_p4est/dg_2d.jl @@ -138,46 +138,46 @@ function prolong2interfaces!(cache, u, # For step = 0, the direction identified by this index is normal to the face. # For step != 0 (1 or -1), the direction identified by this index is tangential to the face. - # Precompute the interpolation indices for the primary element, case: - - # 1) i is the normal direction (constant), j varies along the surface - # => Interpolate in first/normal direction - interp_side_1 = (primary_indices[1] === :begin) ? 1 : 2 - - # 2) j is the normal direction (constant), i varies along the surface - # => Interpolate in second/normal direction - interp_side_2 = (primary_indices[2] === :begin) ? 1 : 2 - # Note that in the current implementation, the interface will be # "aligned at the primary element", i.e., the index of the primary side # will always run forwards. i_primary = i_primary_start j_primary = j_primary_start - for i in eachnode(dg) - for v in eachvariable(equations) - u_primary = zero(eltype(interfaces.u)) - if i_primary_step == 0 - # Figure out if we go forward or backward (i.e., query element orientation) - # to know if we need to use `boundary_interpolation[:, 1]` or `boundary_interpolation[:, 2]`. - # If the first index of the element is `:begin`, this face corresponds to left reference element side (-1). - # If the first index of the element is `:end`, this face corresponds to right reference element side (+1). + + if i_primary_step == 0 + # i is the normal direction (constant), j varies along the surface + # => Interpolate in first/normal direction + interp_side = (primary_indices[1] === :begin) ? 1 : 2 + + for i in eachnode(dg) + for v in eachvariable(equations) + u_primary = zero(eltype(interfaces.u)) for ii in eachnode(dg) u_primary = (u_primary + u[v, ii, j_primary, primary_element] * - boundary_interpolation[ii, interp_side_1]) + boundary_interpolation[ii, interp_side]) end - else # j_primary_step == 0 + interfaces.u[1, v, i, interface] = u_primary + end + j_primary += j_primary_step # incrementing j_primary suffices + end + else # j_primary_step == 0 + # j is the normal direction (constant), i varies along the surface + # => Interpolate in second/normal direction + interp_side = (primary_indices[2] === :begin) ? 1 : 2 + for i in eachnode(dg) + for v in eachvariable(equations) + u_primary = zero(eltype(interfaces.u)) for jj in eachnode(dg) u_primary = (u_primary + u[v, i_primary, jj, primary_element] * - boundary_interpolation[jj, interp_side_2]) + boundary_interpolation[jj, interp_side]) end + interfaces.u[1, v, i, interface] = u_primary end - interfaces.u[1, v, i, interface] = u_primary + i_primary += i_primary_step # incrementing i_primary suffices end - i_primary += i_primary_step - j_primary += j_primary_step end # Interpolate solution data from the secondary element to the interface. @@ -189,37 +189,40 @@ function prolong2interfaces!(cache, u, j_secondary_start, j_secondary_step = index_to_start_step_2d(secondary_indices[2], index_range) - # Precompute the interpolation indices for the secondary element, case: - # 1) i is the normal direction (constant), j varies along the surface - # => Interpolate in first/normal direction - interp_side_1 = (secondary_indices[1] === :begin) ? 1 : 2 - - # 2) j is the normal direction (constant), i varies along the surface - # => Interpolate in second/normal direction - interp_side_2 = (secondary_indices[2] === :begin) ? 1 : 2 - i_secondary = i_secondary_start j_secondary = j_secondary_start - for i in eachnode(dg) - for v in eachvariable(equations) - u_secondary = zero(eltype(interfaces.u)) - if i_secondary_step == 0 + if i_secondary_step == 0 + # i is the normal direction (constant), j varies along the surface + # => Interpolate in first/normal direction + interp_side = (secondary_indices[1] === :begin) ? 1 : 2 + for i in eachnode(dg) + for v in eachvariable(equations) + u_secondary = zero(eltype(interfaces.u)) for ii in eachnode(dg) u_secondary = (u_secondary + u[v, ii, j_secondary, secondary_element] * - boundary_interpolation[ii, interp_side_1]) + boundary_interpolation[ii, interp_side]) end - else # j_secondary_step == 0 + interfaces.u[2, v, i, interface] = u_secondary + end + j_secondary += j_secondary_step # incrementing j_secondary suffices + end + else # j_secondary_step == 0 + # j is the normal direction (constant), i varies along the surface + # => Interpolate in second/normal direction + interp_side = (secondary_indices[2] === :begin) ? 1 : 2 + for i in eachnode(dg) + for v in eachvariable(equations) + u_secondary = zero(eltype(interfaces.u)) for jj in eachnode(dg) u_secondary = (u_secondary + u[v, i_secondary, jj, secondary_element] * - boundary_interpolation[jj, interp_side_2]) + boundary_interpolation[jj, interp_side]) end + interfaces.u[2, v, i, interface] = u_secondary end - interfaces.u[2, v, i, interface] = u_secondary + i_secondary += i_secondary_step # incrementing i_secondary suffices end - i_secondary += i_secondary_step - j_secondary += j_secondary_step end end From 734b221775436d6d9985e6214777d5efb8c7cdf0 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Sun, 15 Mar 2026 09:18:59 +0100 Subject: [PATCH 10/46] Interface normal directions --- src/meshes/p4est_mesh_view.jl | 1 + src/solvers/dgsem_p4est/containers.jl | 69 ++++++++--- src/solvers/dgsem_p4est/containers_2d.jl | 115 ++++++++++++++++++ .../dgsem_p4est/containers_parallel.jl | 4 +- src/solvers/dgsem_p4est/dg_2d.jl | 65 ++++------ src/solvers/dgsem_t8code/containers.jl | 3 +- 6 files changed, 202 insertions(+), 55 deletions(-) diff --git a/src/meshes/p4est_mesh_view.jl b/src/meshes/p4est_mesh_view.jl index 4f0bce4c69b..5e2377dae7a 100644 --- a/src/meshes/p4est_mesh_view.jl +++ b/src/meshes/p4est_mesh_view.jl @@ -97,6 +97,7 @@ function extract_interfaces(mesh::P4estMeshView, interfaces_parent) # Copy relevant entries from parent mesh @views interfaces.u .= interfaces_parent.u[.., mask] + @views interfaces.normal_directions .= interfaces_parent.normal_directions[.., mask] @views interfaces.node_indices .= interfaces_parent.node_indices[.., mask] @views neighbor_ids = interfaces_parent.neighbor_ids[.., mask] diff --git a/src/solvers/dgsem_p4est/containers.jl b/src/solvers/dgsem_p4est/containers.jl index 7ac6febf470..310dc59889c 100644 --- a/src/solvers/dgsem_p4est/containers.jl +++ b/src/solvers/dgsem_p4est/containers.jl @@ -213,22 +213,27 @@ function Adapt.adapt_structure(to, _surface_flux_values) end -mutable struct P4estInterfaceContainer{NDIMS, uEltype <: Real, NDIMSP2, +mutable struct P4estInterfaceContainer{NDIMS, uEltype <: Real, RealT <: Real, + NDIMSP1, NDIMSP2, uArray <: DenseArray{uEltype, NDIMSP2}, + NormalArray <: DenseArray{RealT, NDIMSP1}, IdsMatrix <: DenseMatrix{Int}, IndicesMatrix <: DenseMatrix{NTuple{NDIMS, Symbol}}, uVector <: DenseVector{uEltype}, + NormalVector <: DenseVector{RealT}, IdsVector <: DenseVector{Int}, IndicesVector <: DenseVector{NTuple{NDIMS, Symbol}}} <: AbstractInterfaceContainer - u::uArray # [primary/secondary, variable, i, j, interface] - neighbor_ids::IdsMatrix # [primary/secondary, interface] - node_indices::IndicesMatrix # [primary/secondary, interface] + u::uArray # [primary/secondary, variable, i, j, interface] + normal_directions::NormalArray # [dimension, i, j, interface] + neighbor_ids::IdsMatrix # [primary/secondary, interface] + node_indices::IndicesMatrix # [primary/secondary, interface] # internal `resize!`able storage _u::uVector + _normal_directions::NormalVector _neighbor_ids::IdsVector _node_indices::IndicesVector end @@ -244,7 +249,7 @@ end # See explanation of Base.resize! for the element container function Base.resize!(interfaces::P4estInterfaceContainer, capacity) - @unpack _u, _neighbor_ids, _node_indices = interfaces + @unpack _u, _normal_directions, _neighbor_ids, _node_indices = interfaces n_dims = ndims(interfaces) n_nodes = size(interfaces.u, 3) @@ -256,6 +261,12 @@ function Base.resize!(interfaces::P4estInterfaceContainer, capacity) (2, n_variables, ntuple(_ -> n_nodes, n_dims - 1)..., capacity)) + resize!(_normal_directions, n_dims * n_nodes^(n_dims - 1) * capacity) + interfaces.normal_directions = unsafe_wrap(ArrayType, pointer(_normal_directions), + (n_dims, + ntuple(_ -> n_nodes, n_dims - 1)..., + capacity)) + resize!(_neighbor_ids, 2 * capacity) interfaces.neighbor_ids = unsafe_wrap(ArrayType, pointer(_neighbor_ids), (2, capacity)) @@ -271,6 +282,7 @@ end function init_interfaces(mesh::Union{P4estMesh, P4estMeshView, T8codeMesh}, equations, basis, elements) NDIMS = ndims(elements) + RealT = real(mesh) uEltype = eltype(elements) # Initialize container @@ -283,34 +295,49 @@ function init_interfaces(mesh::Union{P4estMesh, P4estMeshView, T8codeMesh}, equa (2, nvariables(equations), ntuple(_ -> nnodes(basis), NDIMS - 1)..., n_interfaces)) + _normal_directions = Vector{RealT}(undef, + NDIMS * nnodes(basis)^(NDIMS - 1) * + n_interfaces) + normal_directions = unsafe_wrap(Array, pointer(_normal_directions), + (NDIMS, + ntuple(_ -> nnodes(basis), NDIMS - 1)..., + n_interfaces)) + _neighbor_ids = Vector{Int}(undef, 2 * n_interfaces) neighbor_ids = unsafe_wrap(Array, pointer(_neighbor_ids), (2, n_interfaces)) _node_indices = Vector{NTuple{NDIMS, Symbol}}(undef, 2 * n_interfaces) node_indices = unsafe_wrap(Array, pointer(_node_indices), (2, n_interfaces)) - interfaces = P4estInterfaceContainer{NDIMS, uEltype, NDIMS + 2, - typeof(u), typeof(neighbor_ids), - typeof(node_indices), typeof(_u), + interfaces = P4estInterfaceContainer{NDIMS, uEltype, RealT, + NDIMS + 1, NDIMS + 2, + typeof(u), typeof(normal_directions), + typeof(neighbor_ids), typeof(node_indices), + typeof(_u), typeof(_normal_directions), typeof(_neighbor_ids), typeof(_node_indices)}(u, + normal_directions, neighbor_ids, node_indices, _u, + _normal_directions, _neighbor_ids, _node_indices) - init_interfaces!(interfaces, mesh) + init_interfaces!(interfaces, elements, mesh, basis) return interfaces end -function init_interfaces!(interfaces, mesh::Union{P4estMesh, P4estMeshView}) +function init_interfaces!(interfaces, elements, + mesh::Union{P4estMesh, P4estMeshView}, basis) init_surfaces!(interfaces, nothing, nothing, mesh) + init_normal_directions!(interfaces, basis, elements) return interfaces end function Adapt.parent_type(::Type{<:P4estInterfaceContainer{<:Any, <:Any, <:Any, + <:Any, <:Any, ArrayT}}) where {ArrayT} return ArrayT end @@ -319,10 +346,13 @@ end function Adapt.adapt_structure(to, interfaces::P4estInterfaceContainer) # Adapt underlying storage _u = adapt(to, interfaces._u) + _normal_directions = adapt(to, interfaces._normal_directions) _neighbor_ids = adapt(to, interfaces._neighbor_ids) _node_indices = adapt(to, interfaces._node_indices) # Wrap arrays again u = unsafe_wrap_or_alloc(to, _u, size(interfaces.u)) + normal_directions = unsafe_wrap_or_alloc(to, _normal_directions, + size(interfaces.normal_directions)) neighbor_ids = unsafe_wrap_or_alloc(to, _neighbor_ids, size(interfaces.neighbor_ids)) node_indices = unsafe_wrap_or_alloc(to, _node_indices, @@ -331,11 +361,16 @@ function Adapt.adapt_structure(to, interfaces::P4estInterfaceContainer) NDIMS = ndims(interfaces) new_type_params = (NDIMS, eltype(_u), - NDIMS + 2, - typeof(u), typeof(neighbor_ids), typeof(node_indices), - typeof(_u), typeof(_neighbor_ids), typeof(_node_indices)) - return P4estInterfaceContainer{new_type_params...}(u, neighbor_ids, node_indices, - _u, _neighbor_ids, _node_indices) + eltype(_normal_directions), + NDIMS + 1, NDIMS + 2, + typeof(u), typeof(normal_directions), + typeof(neighbor_ids), typeof(node_indices), + typeof(_u), typeof(_normal_directions), + typeof(_neighbor_ids), typeof(_node_indices)) + return P4estInterfaceContainer{new_type_params...}(u, normal_directions, + neighbor_ids, node_indices, + _u, _normal_directions, + _neighbor_ids, _node_indices) end mutable struct P4estBoundaryContainer{NDIMS, uEltype <: Real, NDIMSP1, @@ -644,6 +679,10 @@ function reinitialize_containers!(mesh::P4estMesh, equations, dg::DGSEM, cache) # resize interfaces container @unpack interfaces = cache resize!(interfaces, required.interfaces) + # Call `init_normal_directions!` only since `init_interfaces!` calls also + # `init_surfaces!` which iterates over the mesh, which should be done only once + # after boundaries and mortars are resized. + init_normal_directions!(interfaces, dg.basis, elements) # resize boundaries container @unpack boundaries = cache diff --git a/src/solvers/dgsem_p4est/containers_2d.jl b/src/solvers/dgsem_p4est/containers_2d.jl index 9f6d31f7ed0..6cfdb6e7eb5 100644 --- a/src/solvers/dgsem_p4est/containers_2d.jl +++ b/src/solvers/dgsem_p4est/containers_2d.jl @@ -87,6 +87,121 @@ function calc_node_coordinates!(node_coordinates, return node_coordinates end +# For Gauss-Lobatto-Legendre (GLL) nodes, the interface normals are +# simply taken from the element face nodes. +function init_normal_directions!(interfaces::P4estInterfaceContainer{2}, + basis::LobattoLegendreBasis, + elements::P4estElementContainer{2}) + @unpack neighbor_ids, node_indices, normal_directions = interfaces + @unpack contravariant_vectors = elements + index_range = eachnode(basis) + index_begin = first(index_range) + index_end = last(index_range) + + for interface in axes(neighbor_ids, 2) + primary_element = neighbor_ids[1, interface] + primary_indices = node_indices[1, interface] + primary_direction = indices2direction(primary_indices) + + i_primary, i_primary_step = index_to_start_step_2d(primary_indices[1], + index_begin, + index_end) + j_primary, j_primary_step = index_to_start_step_2d(primary_indices[2], + index_begin, + index_end) + + for i in index_range + # Retrieve normal direction at element/volume nodes + normal_direction = get_normal_direction(primary_direction, + contravariant_vectors, + i_primary, j_primary, + primary_element) + normal_directions[1, i, interface] = normal_direction[1] + normal_directions[2, i, interface] = normal_direction[2] + + i_primary += i_primary_step + j_primary += j_primary_step + end + end + + return interfaces +end + +function init_normal_directions!(interfaces::P4estInterfaceContainer{2}, + basis::GaussLegendreBasis, + elements::P4estElementContainer{2}) + @unpack neighbor_ids, node_indices, normal_directions = interfaces + @unpack contravariant_vectors = elements + @unpack boundary_interpolation = basis + index_range = eachnode(basis) + index_begin = first(index_range) + index_end = last(index_range) + + for interface in axes(neighbor_ids, 2) + primary_element = neighbor_ids[1, interface] + primary_indices = node_indices[1, interface] + primary_direction = indices2direction(primary_indices) + + i_primary_start, i_primary_step = index_to_start_step_2d(primary_indices[1], + index_begin, + index_end) + j_primary_start, j_primary_step = index_to_start_step_2d(primary_indices[2], + index_begin, + index_end) + + i_primary = i_primary_start + j_primary = j_primary_start + + if i_primary_step == 0 + # i is the normal direction (constant), j varies along the surface + # => Interpolate in first/normal direction + # Interpolation side is governed by element orientation + side = interpolation_side(primary_indices[1]) + for i in index_range + normal_1 = zero(eltype(normal_directions)) + normal_2 = zero(eltype(normal_directions)) + for ii in index_range + factor = boundary_interpolation[ii, side] + # Retrieve normal directions at element/volume nodes + normal_direction = get_normal_direction(primary_direction, + contravariant_vectors, + ii, j_primary, + primary_element) + normal_1 = normal_1 + normal_direction[1] * factor + normal_2 = normal_2 + normal_direction[2] * factor + end + normal_directions[1, i, interface] = normal_1 + normal_directions[2, i, interface] = normal_2 + j_primary += j_primary_step # incrementing j_primary suffices (i_primary_step = 0) + end + else # j_primary_step == 0 + # j is the normal direction (constant), i varies along the surface + # => Interpolate in second/normal direction + # Interpolation side is governed by element orientation + side = interpolation_side(primary_indices[2]) + for i in index_range + normal_1 = zero(eltype(normal_directions)) + normal_2 = zero(eltype(normal_directions)) + for jj in index_range + factor = boundary_interpolation[jj, side] + # Retrieve normal directions at element/volume nodes + normal_direction = get_normal_direction(primary_direction, + contravariant_vectors, + i_primary, jj, + primary_element) + normal_1 = normal_1 + normal_direction[1] * factor + normal_2 = normal_2 + normal_direction[2] * factor + end + normal_directions[1, i, interface] = normal_1 + normal_directions[2, i, interface] = normal_2 + i_primary += i_primary_step # incrementing i_primary suffices (j_primary_step = 0) + end + end + end + + return interfaces +end + # Initialize node_indices of interface container @inline function init_interface_node_indices!(interfaces::P4estInterfaceContainer{2}, faces, orientation, interface_id) diff --git a/src/solvers/dgsem_p4est/containers_parallel.jl b/src/solvers/dgsem_p4est/containers_parallel.jl index 82cf2bcbf56..ea29fb3cdae 100644 --- a/src/solvers/dgsem_p4est/containers_parallel.jl +++ b/src/solvers/dgsem_p4est/containers_parallel.jl @@ -252,8 +252,10 @@ end # Overload init! function for regular interfaces, regular mortars and boundaries since they must # call the appropriate init_surfaces! function for parallel p4est meshes -function init_interfaces!(interfaces, mesh::P4estMeshParallel) +function init_interfaces!(interfaces, elements, + mesh::P4estMeshParallel, basis) init_surfaces!(interfaces, nothing, nothing, nothing, nothing, mesh) + init_normal_directions!(interfaces, basis, elements) return interfaces end diff --git a/src/solvers/dgsem_p4est/dg_2d.jl b/src/solvers/dgsem_p4est/dg_2d.jl index 61a865c4f37..42f3c6b8c2f 100644 --- a/src/solvers/dgsem_p4est/dg_2d.jl +++ b/src/solvers/dgsem_p4est/dg_2d.jl @@ -26,6 +26,7 @@ function create_cache(mesh::Union{P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}} return cache end +# index_to_start_step_2d(index::Symbol, index_begin, index_end) # index_to_start_step_2d(index::Symbol, index_range) # # Given a symbolic `index` and an `indexrange` (usually `eachnode(dg)`), @@ -47,10 +48,7 @@ end # i_volume += i_volume_step # j_volume += j_volume_step # end -@inline function index_to_start_step_2d(index::Symbol, index_range) - index_begin = first(index_range) - index_end = last(index_range) - +@inline function index_to_start_step_2d(index::Symbol, index_begin, index_end) if index === :begin return index_begin, 0 elseif index === :end @@ -61,6 +59,18 @@ end return index_end, -1 end end +@inline function index_to_start_step_2d(index::Symbol, index_range) + index_begin = first(index_range) + index_end = last(index_range) + + return index_to_start_step_2d(index, index_begin, index_end) +end + +# Infer interpolation side, i.e., left (1) or right (2) for an element. +# Required for boundary interpolation with Gauss-Legendre nodes. +@inline function interpolation_side(index) + return (index === :begin) ? 1 : 2 +end function prolong2interfaces!(cache, u, mesh::Union{P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}}, @@ -148,15 +158,15 @@ function prolong2interfaces!(cache, u, if i_primary_step == 0 # i is the normal direction (constant), j varies along the surface # => Interpolate in first/normal direction - interp_side = (primary_indices[1] === :begin) ? 1 : 2 - + # Interpolation side is governed by element orientation + side = interpolation_side(primary_indices[1]) for i in eachnode(dg) for v in eachvariable(equations) u_primary = zero(eltype(interfaces.u)) for ii in eachnode(dg) u_primary = (u_primary + u[v, ii, j_primary, primary_element] * - boundary_interpolation[ii, interp_side]) + boundary_interpolation[ii, side]) end interfaces.u[1, v, i, interface] = u_primary end @@ -165,14 +175,15 @@ function prolong2interfaces!(cache, u, else # j_primary_step == 0 # j is the normal direction (constant), i varies along the surface # => Interpolate in second/normal direction - interp_side = (primary_indices[2] === :begin) ? 1 : 2 + # Interpolation side is governed by element orientation + side = interpolation_side(primary_indices[2]) for i in eachnode(dg) for v in eachvariable(equations) u_primary = zero(eltype(interfaces.u)) for jj in eachnode(dg) u_primary = (u_primary + u[v, i_primary, jj, primary_element] * - boundary_interpolation[jj, interp_side]) + boundary_interpolation[jj, side]) end interfaces.u[1, v, i, interface] = u_primary end @@ -192,32 +203,28 @@ function prolong2interfaces!(cache, u, i_secondary = i_secondary_start j_secondary = j_secondary_start if i_secondary_step == 0 - # i is the normal direction (constant), j varies along the surface - # => Interpolate in first/normal direction - interp_side = (secondary_indices[1] === :begin) ? 1 : 2 + side = interpolation_side(secondary_indices[1]) for i in eachnode(dg) for v in eachvariable(equations) u_secondary = zero(eltype(interfaces.u)) for ii in eachnode(dg) u_secondary = (u_secondary + u[v, ii, j_secondary, secondary_element] * - boundary_interpolation[ii, interp_side]) + boundary_interpolation[ii, side]) end interfaces.u[2, v, i, interface] = u_secondary end j_secondary += j_secondary_step # incrementing j_secondary suffices end else # j_secondary_step == 0 - # j is the normal direction (constant), i varies along the surface - # => Interpolate in second/normal direction - interp_side = (secondary_indices[2] === :begin) ? 1 : 2 + side = interpolation_side(secondary_indices[2]) for i in eachnode(dg) for v in eachvariable(equations) u_secondary = zero(eltype(interfaces.u)) for jj in eachnode(dg) u_secondary = (u_secondary + u[v, i_secondary, jj, secondary_element] * - boundary_interpolation[jj, interp_side]) + boundary_interpolation[jj, side]) end interfaces.u[2, v, i, interface] = u_secondary end @@ -234,8 +241,7 @@ function calc_interface_flux!(surface_flux_values, T8codeMesh{2}}, have_nonconservative_terms, equations, surface_integral, dg::DG, cache) - @unpack neighbor_ids, node_indices = cache.interfaces - @unpack contravariant_vectors = cache.elements + @unpack neighbor_ids, node_indices, normal_directions = cache.interfaces index_range = eachnode(dg) index_end = last(index_range) @@ -245,15 +251,6 @@ function calc_interface_flux!(surface_flux_values, primary_indices = node_indices[1, interface] primary_direction = indices2direction(primary_indices) - # Create the local i,j indexing on the primary element used to pull normal direction information - i_primary_start, i_primary_step = index_to_start_step_2d(primary_indices[1], - index_range) - j_primary_start, j_primary_step = index_to_start_step_2d(primary_indices[2], - index_range) - - i_primary = i_primary_start - j_primary = j_primary_start - # Get element and side index information on the secondary element secondary_element = neighbor_ids[2, interface] secondary_indices = node_indices[2, interface] @@ -271,13 +268,8 @@ function calc_interface_flux!(surface_flux_values, end for node in eachnode(dg) - # Get the normal direction on the primary element. - # Contravariant vectors at interfaces in negative coordinate direction - # are pointing inwards. This is handled by `get_normal_direction`. - normal_direction = get_normal_direction(primary_direction, - contravariant_vectors, - i_primary, j_primary, - primary_element) + normal_direction = SVector(normal_directions[1, node, interface], + normal_directions[2, node, interface]) calc_interface_flux!(surface_flux_values, mesh, have_nonconservative_terms, equations, @@ -286,9 +278,6 @@ function calc_interface_flux!(surface_flux_values, node, primary_direction, primary_element, node_secondary, secondary_direction, secondary_element) - # Increment primary element indices to pull the normal direction - i_primary += i_primary_step - j_primary += j_primary_step # Increment the surface node index along the secondary element node_secondary += node_secondary_step end diff --git a/src/solvers/dgsem_t8code/containers.jl b/src/solvers/dgsem_t8code/containers.jl index dea0bf1a9f1..56db011f3b1 100644 --- a/src/solvers/dgsem_t8code/containers.jl +++ b/src/solvers/dgsem_t8code/containers.jl @@ -54,7 +54,8 @@ function count_required_surfaces(mesh::T8codeMesh) end # Compatibility to `dgsem_p4est/containers.jl`. -function init_interfaces!(interfaces, mesh::T8codeMesh) +function init_interfaces!(interfaces, elements, + mesh::T8codeMesh, basis) # Already computed. Do nothing. return nothing end From f17dcd31bf08344767010cfdc27b2f76e4f52fc2 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Sun, 15 Mar 2026 09:23:33 +0100 Subject: [PATCH 11/46] example --- .../elixir_advection_flag_gauss_legendre.jl | 49 +++++++++++++++++++ src/solvers/dgsem_p4est/dg_2d.jl | 1 + test/test_p4est_2d.jl | 9 ++++ 3 files changed, 59 insertions(+) create mode 100644 examples/p4est_2d_dgsem/elixir_advection_flag_gauss_legendre.jl diff --git a/examples/p4est_2d_dgsem/elixir_advection_flag_gauss_legendre.jl b/examples/p4est_2d_dgsem/elixir_advection_flag_gauss_legendre.jl new file mode 100644 index 00000000000..cbe37f087e3 --- /dev/null +++ b/examples/p4est_2d_dgsem/elixir_advection_flag_gauss_legendre.jl @@ -0,0 +1,49 @@ +using OrdinaryDiffEqLowStorageRK +using Trixi + +############################################################################### +# semidiscretization of the linear advection equation + +advection_velocity = (0.2, -0.7) +equations = LinearScalarAdvectionEquation2D(advection_velocity) + +solver = DGSEM(polydeg = 4, surface_flux = flux_lax_friedrichs, + basis_type = GaussLegendreBasis) + +# Deformed rectangle that looks like a waving flag, +# lower and upper faces are sinus curves, left and right are vertical lines. +f1(s) = SVector(-1.0, s - 1.0) +f2(s) = SVector(1.0, s + 1.0) +f3(s) = SVector(s, -1.0 + sin(0.5 * pi * s)) +f4(s) = SVector(s, 1.0 + sin(0.5 * pi * s)) + +# Create P4estMesh with 3 x 2 trees and 6 x 4 elements, +# approximate the geometry with a smaller polydeg for testing. +trees_per_dimension = (3, 2) +mesh = P4estMesh(trees_per_dimension, polydeg = 3, + faces = (f1, f2, f3, f4), + periodicity = true, + initial_refinement_level = 1) + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_convergence_test, + solver; + boundary_conditions = boundary_condition_periodic) + +############################################################################### +# ODE solvers, callbacks etc. + +ode = semidiscretize(semi, (0.0, 0.2)) + +summary_callback = SummaryCallback() +analysis_callback = AnalysisCallback(semi, interval = 100) +stepsize_callback = StepsizeCallback(cfl = 0.8) + +callbacks = CallbackSet(summary_callback, #analysis_callback, + stepsize_callback) + +############################################################################### +# run the simulation + +sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false); + dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + ode_default_options()..., callback = callbacks); diff --git a/src/solvers/dgsem_p4est/dg_2d.jl b/src/solvers/dgsem_p4est/dg_2d.jl index 42f3c6b8c2f..544bcff4eef 100644 --- a/src/solvers/dgsem_p4est/dg_2d.jl +++ b/src/solvers/dgsem_p4est/dg_2d.jl @@ -268,6 +268,7 @@ function calc_interface_flux!(surface_flux_values, end for node in eachnode(dg) + # This seems to be faster than `@views normal_direction = normal_directions[:, node, interface]` normal_direction = SVector(normal_directions[1, node, interface], normal_directions[2, node, interface]) diff --git a/test/test_p4est_2d.jl b/test/test_p4est_2d.jl index b030dfe4fe4..41e6a13bb61 100644 --- a/test/test_p4est_2d.jl +++ b/test/test_p4est_2d.jl @@ -73,6 +73,15 @@ end @test_allocations(Trixi.rhs!, semi, sol, 1000) end +@trixi_testset "elixir_advection_flag_gauss_legendre.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_advection_flag_gauss_legendre.jl"), + l2=[0.00047332539212915036], linf=[0.0021788728950808967]) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + @test_allocations(Trixi.rhs!, semi, sol, 1000) +end + @trixi_testset "elixir_advection_unstructured_flag.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_unstructured_flag.jl"), l2=[0.0005379687442422346], From d69e57061955c06acd5822fa85d4f8041e2eb1ea Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Sun, 15 Mar 2026 09:29:21 +0100 Subject: [PATCH 12/46] comment --- src/solvers/dgsem_p4est/containers_2d.jl | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/solvers/dgsem_p4est/containers_2d.jl b/src/solvers/dgsem_p4est/containers_2d.jl index 6cfdb6e7eb5..a2b1e967bb5 100644 --- a/src/solvers/dgsem_p4est/containers_2d.jl +++ b/src/solvers/dgsem_p4est/containers_2d.jl @@ -127,6 +127,8 @@ function init_normal_directions!(interfaces::P4estInterfaceContainer{2}, return interfaces end +# For Gauss-Legendre (GL) nodes, the interface normals are +# computed from interpolation of the volume node normals to the surface nodes. function init_normal_directions!(interfaces::P4estInterfaceContainer{2}, basis::GaussLegendreBasis, elements::P4estElementContainer{2}) From 63eb6966178401e2f12972d38b71d6a06dc8bcf3 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Sun, 15 Mar 2026 09:48:12 +0100 Subject: [PATCH 13/46] 3d dummy version --- src/solvers/dgsem_p4est/containers_3d.jl | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/src/solvers/dgsem_p4est/containers_3d.jl b/src/solvers/dgsem_p4est/containers_3d.jl index 9d343c76f02..c01e04d6f5e 100644 --- a/src/solvers/dgsem_p4est/containers_3d.jl +++ b/src/solvers/dgsem_p4est/containers_3d.jl @@ -76,6 +76,13 @@ function calc_node_coordinates!(node_coordinates, return node_coordinates end +# Not yet implmented and needed for 3D +function init_normal_directions!(interfaces::P4estInterfaceContainer{3}, + basis::LobattoLegendreBasis, + elements::P4estElementContainer{3}) + return nothing +end + # Initialize node_indices of interface container @inline function init_interface_node_indices!(interfaces::P4estInterfaceContainer{3}, faces, orientation, interface_id) From 4672a79cee48717b603ccfcf41a5cd7514d3f0de Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Sun, 15 Mar 2026 09:49:21 +0100 Subject: [PATCH 14/46] typo --- src/solvers/dgsem_p4est/containers_3d.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/solvers/dgsem_p4est/containers_3d.jl b/src/solvers/dgsem_p4est/containers_3d.jl index c01e04d6f5e..55972861e64 100644 --- a/src/solvers/dgsem_p4est/containers_3d.jl +++ b/src/solvers/dgsem_p4est/containers_3d.jl @@ -76,7 +76,7 @@ function calc_node_coordinates!(node_coordinates, return node_coordinates end -# Not yet implmented and needed for 3D +# Not yet implemented and needed for 3D function init_normal_directions!(interfaces::P4estInterfaceContainer{3}, basis::LobattoLegendreBasis, elements::P4estElementContainer{3}) From 920037b46eef5d4c026f6b0997361320e6ee2521 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Sun, 15 Mar 2026 10:10:40 +0100 Subject: [PATCH 15/46] fix --- src/solvers/dgsem_p4est/dg_2d.jl | 3 +- src/solvers/dgsem_t8code/dg.jl | 2 + src/solvers/dgsem_t8code/dg_2d.jl | 72 +++++++++++++++++++++++++++++++ 3 files changed, 75 insertions(+), 2 deletions(-) create mode 100644 src/solvers/dgsem_t8code/dg_2d.jl diff --git a/src/solvers/dgsem_p4est/dg_2d.jl b/src/solvers/dgsem_p4est/dg_2d.jl index 544bcff4eef..b6be0e7a947 100644 --- a/src/solvers/dgsem_p4est/dg_2d.jl +++ b/src/solvers/dgsem_p4est/dg_2d.jl @@ -237,8 +237,7 @@ function prolong2interfaces!(cache, u, end function calc_interface_flux!(surface_flux_values, - mesh::Union{P4estMesh{2}, P4estMeshView{2}, - T8codeMesh{2}}, + mesh::Union{P4estMesh{2}, P4estMeshView{2}}, have_nonconservative_terms, equations, surface_integral, dg::DG, cache) @unpack neighbor_ids, node_indices, normal_directions = cache.interfaces diff --git a/src/solvers/dgsem_t8code/dg.jl b/src/solvers/dgsem_t8code/dg.jl index ce35d524f8b..acab6330097 100644 --- a/src/solvers/dgsem_t8code/dg.jl +++ b/src/solvers/dgsem_t8code/dg.jl @@ -33,5 +33,7 @@ include("containers_2d.jl") include("containers_3d.jl") include("containers_parallel.jl") + +include("dg_2d.jl") include("dg_parallel.jl") end # @muladd diff --git a/src/solvers/dgsem_t8code/dg_2d.jl b/src/solvers/dgsem_t8code/dg_2d.jl new file mode 100644 index 00000000000..308b268f196 --- /dev/null +++ b/src/solvers/dgsem_t8code/dg_2d.jl @@ -0,0 +1,72 @@ +@muladd begin +#! format: noindent + +# Use for `T8CodeMesh` the Gauss-Lobatto-Legendre (GLL) hard-coded variant, +# i.e., take interface normals from the outer volume nodes. +function calc_interface_flux!(surface_flux_values, + mesh::T8codeMesh{2}, + have_nonconservative_terms, + equations, surface_integral, dg::DG, cache) + @unpack neighbor_ids, node_indices = cache.interfaces + @unpack contravariant_vectors = cache.elements + index_range = eachnode(dg) + index_end = last(index_range) + + @threaded for interface in eachinterface(dg, cache) + # Get element and side index information on the primary element + primary_element = neighbor_ids[1, interface] + primary_indices = node_indices[1, interface] + primary_direction = indices2direction(primary_indices) + + # Create the local i,j indexing on the primary element used to pull normal direction information + i_primary_start, i_primary_step = index_to_start_step_2d(primary_indices[1], + index_range) + j_primary_start, j_primary_step = index_to_start_step_2d(primary_indices[2], + index_range) + + i_primary = i_primary_start + j_primary = j_primary_start + + # Get element and side index information on the secondary element + secondary_element = neighbor_ids[2, interface] + secondary_indices = node_indices[2, interface] + secondary_direction = indices2direction(secondary_indices) + + # Initiate the secondary index to be used in the surface for loop. + # This index on the primary side will always run forward but + # the secondary index might need to run backwards for flipped sides. + if :i_backward in secondary_indices + node_secondary = index_end + node_secondary_step = -1 + else + node_secondary = 1 + node_secondary_step = 1 + end + + for node in eachnode(dg) + # Get the normal direction on the primary element. + # Contravariant vectors at interfaces in negative coordinate direction + # are pointing inwards. This is handled by `get_normal_direction`. + normal_direction = get_normal_direction(primary_direction, + contravariant_vectors, + i_primary, j_primary, + primary_element) + + calc_interface_flux!(surface_flux_values, mesh, have_nonconservative_terms, + equations, + surface_integral, dg, cache, + interface, normal_direction, + node, primary_direction, primary_element, + node_secondary, secondary_direction, secondary_element) + + # Increment primary element indices to pull the normal direction + i_primary += i_primary_step + j_primary += j_primary_step + # Increment the surface node index along the secondary element + node_secondary += node_secondary_step + end + end + + return nothing +end +end # @muladd From 2b797a4c40a62b9d0a3f7e691c7b55f634395e6b Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Sun, 15 Mar 2026 10:14:06 +0100 Subject: [PATCH 16/46] comment --- src/solvers/dgsem_p4est/containers_2d.jl | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/solvers/dgsem_p4est/containers_2d.jl b/src/solvers/dgsem_p4est/containers_2d.jl index a2b1e967bb5..1934b44e84a 100644 --- a/src/solvers/dgsem_p4est/containers_2d.jl +++ b/src/solvers/dgsem_p4est/containers_2d.jl @@ -151,6 +151,10 @@ function init_normal_directions!(interfaces::P4estInterfaceContainer{2}, index_begin, index_end) + # The index direction is identified based on `{i,j}_{primary, secondary}_step`. + # For step = 0, the direction identified by this index is normal to the face. + # For step != 0 (1 or -1), the direction identified by this index is tangential to the face. + i_primary = i_primary_start j_primary = j_primary_start From 9cce03b3c783efa6bcb3314a8c42d0094c97f121 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Sun, 15 Mar 2026 10:19:33 +0100 Subject: [PATCH 17/46] rev --- src/solvers/dgsem_tree/dg_1d.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/solvers/dgsem_tree/dg_1d.jl b/src/solvers/dgsem_tree/dg_1d.jl index 7f98883a2b5..6bc81900953 100644 --- a/src/solvers/dgsem_tree/dg_1d.jl +++ b/src/solvers/dgsem_tree/dg_1d.jl @@ -701,7 +701,7 @@ function calc_surface_integral!(du, u, mesh::Union{TreeMesh{1}, StructuredMesh{1 return nothing end -function calc_surface_integral!(du, u, mesh::TreeMesh{1}, +function calc_surface_integral!(du, u, mesh::Union{TreeMesh{1}, StructuredMesh{1}}, equations, surface_integral::SurfaceIntegralWeakForm, dg::DGSEM{<:GaussLegendreBasis}, cache) @unpack boundary_interpolation_inverse_weights = dg.basis From 6ba5bc755661d0a13de2f9cadaef47e63175f994 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Sun, 15 Mar 2026 12:49:31 +0100 Subject: [PATCH 18/46] db --- src/solvers/dgsem_p4est/containers.jl | 11 ++++++----- src/solvers/dgsem_p4est/containers_parallel.jl | 4 +++- 2 files changed, 9 insertions(+), 6 deletions(-) diff --git a/src/solvers/dgsem_p4est/containers.jl b/src/solvers/dgsem_p4est/containers.jl index 310dc59889c..9dea096b055 100644 --- a/src/solvers/dgsem_p4est/containers.jl +++ b/src/solvers/dgsem_p4est/containers.jl @@ -679,10 +679,6 @@ function reinitialize_containers!(mesh::P4estMesh, equations, dg::DGSEM, cache) # resize interfaces container @unpack interfaces = cache resize!(interfaces, required.interfaces) - # Call `init_normal_directions!` only since `init_interfaces!` calls also - # `init_surfaces!` which iterates over the mesh, which should be done only once - # after boundaries and mortars are resized. - init_normal_directions!(interfaces, dg.basis, elements) # resize boundaries container @unpack boundaries = cache @@ -694,7 +690,12 @@ function reinitialize_containers!(mesh::P4estMesh, equations, dg::DGSEM, cache) # re-initialize containers together to reduce # the number of iterations over the mesh in `p4est` - return init_surfaces!(interfaces, mortars, boundaries, mesh) + init_surfaces!(interfaces, mortars, boundaries, mesh) + + # Normal directions require that `node_indices` have been initialized + init_normal_directions!(interfaces, dg.basis, elements) + + return nothing end # A helper struct used in initialization methods below diff --git a/src/solvers/dgsem_p4est/containers_parallel.jl b/src/solvers/dgsem_p4est/containers_parallel.jl index ea29fb3cdae..e67ad26606d 100644 --- a/src/solvers/dgsem_p4est/containers_parallel.jl +++ b/src/solvers/dgsem_p4est/containers_parallel.jl @@ -322,7 +322,9 @@ function reinitialize_containers!(mesh::P4estMeshParallel, equations, dg::DGSEM, # re-initialize and distribute normal directions of MPI mortars; requires MPI communication, so # the MPI cache must be re-initialized before init_normal_directions!(mpi_mortars, dg.basis, elements) - return exchange_normal_directions!(mpi_mortars, mpi_cache, mesh, nnodes(dg)) + exchange_normal_directions!(mpi_mortars, mpi_cache, mesh, nnodes(dg)) + + return nothing end # A helper struct used in initialization methods below From 4a36ed973d3de51ce17d6e1c3bf4aff07dc2ee84 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Sun, 15 Mar 2026 12:55:59 +0100 Subject: [PATCH 19/46] rev --- src/solvers/dgsem_p4est/dg_2d.jl | 3 +- src/solvers/dgsem_t8code/containers.jl | 3 + .../dgsem_t8code/containers_parallel.jl | 3 + src/solvers/dgsem_t8code/dg.jl | 1 - src/solvers/dgsem_t8code/dg_2d.jl | 72 ------------------- 5 files changed, 8 insertions(+), 74 deletions(-) delete mode 100644 src/solvers/dgsem_t8code/dg_2d.jl diff --git a/src/solvers/dgsem_p4est/dg_2d.jl b/src/solvers/dgsem_p4est/dg_2d.jl index b6be0e7a947..544bcff4eef 100644 --- a/src/solvers/dgsem_p4est/dg_2d.jl +++ b/src/solvers/dgsem_p4est/dg_2d.jl @@ -237,7 +237,8 @@ function prolong2interfaces!(cache, u, end function calc_interface_flux!(surface_flux_values, - mesh::Union{P4estMesh{2}, P4estMeshView{2}}, + mesh::Union{P4estMesh{2}, P4estMeshView{2}, + T8codeMesh{2}}, have_nonconservative_terms, equations, surface_integral, dg::DG, cache) @unpack neighbor_ids, node_indices, normal_directions = cache.interfaces diff --git a/src/solvers/dgsem_t8code/containers.jl b/src/solvers/dgsem_t8code/containers.jl index 56db011f3b1..ec5c6c9e9dd 100644 --- a/src/solvers/dgsem_t8code/containers.jl +++ b/src/solvers/dgsem_t8code/containers.jl @@ -28,6 +28,9 @@ function reinitialize_containers!(mesh::T8codeMesh, equations, dg::DGSEM, cache) fill_mesh_info!(mesh, interfaces, mortars, boundaries, mesh.boundary_names) + # Normal directions require that `node_indices` have been initialized + init_normal_directions!(interfaces, dg.basis, elements) + return nothing end diff --git a/src/solvers/dgsem_t8code/containers_parallel.jl b/src/solvers/dgsem_t8code/containers_parallel.jl index 65d4209a2c6..27586c5b841 100644 --- a/src/solvers/dgsem_t8code/containers_parallel.jl +++ b/src/solvers/dgsem_t8code/containers_parallel.jl @@ -52,6 +52,9 @@ function reinitialize_containers!(mesh::T8codeMeshParallel, equations, dg::DGSEM init_normal_directions!(mpi_mortars, dg.basis, elements) exchange_normal_directions!(mpi_mortars, mpi_cache, mesh, nnodes(dg)) + # Normal directions require that `node_indices` have been initialized + init_normal_directions!(interfaces, dg.basis, elements) + return nothing end diff --git a/src/solvers/dgsem_t8code/dg.jl b/src/solvers/dgsem_t8code/dg.jl index acab6330097..47169c5d8e7 100644 --- a/src/solvers/dgsem_t8code/dg.jl +++ b/src/solvers/dgsem_t8code/dg.jl @@ -34,6 +34,5 @@ include("containers_3d.jl") include("containers_parallel.jl") -include("dg_2d.jl") include("dg_parallel.jl") end # @muladd diff --git a/src/solvers/dgsem_t8code/dg_2d.jl b/src/solvers/dgsem_t8code/dg_2d.jl deleted file mode 100644 index 308b268f196..00000000000 --- a/src/solvers/dgsem_t8code/dg_2d.jl +++ /dev/null @@ -1,72 +0,0 @@ -@muladd begin -#! format: noindent - -# Use for `T8CodeMesh` the Gauss-Lobatto-Legendre (GLL) hard-coded variant, -# i.e., take interface normals from the outer volume nodes. -function calc_interface_flux!(surface_flux_values, - mesh::T8codeMesh{2}, - have_nonconservative_terms, - equations, surface_integral, dg::DG, cache) - @unpack neighbor_ids, node_indices = cache.interfaces - @unpack contravariant_vectors = cache.elements - index_range = eachnode(dg) - index_end = last(index_range) - - @threaded for interface in eachinterface(dg, cache) - # Get element and side index information on the primary element - primary_element = neighbor_ids[1, interface] - primary_indices = node_indices[1, interface] - primary_direction = indices2direction(primary_indices) - - # Create the local i,j indexing on the primary element used to pull normal direction information - i_primary_start, i_primary_step = index_to_start_step_2d(primary_indices[1], - index_range) - j_primary_start, j_primary_step = index_to_start_step_2d(primary_indices[2], - index_range) - - i_primary = i_primary_start - j_primary = j_primary_start - - # Get element and side index information on the secondary element - secondary_element = neighbor_ids[2, interface] - secondary_indices = node_indices[2, interface] - secondary_direction = indices2direction(secondary_indices) - - # Initiate the secondary index to be used in the surface for loop. - # This index on the primary side will always run forward but - # the secondary index might need to run backwards for flipped sides. - if :i_backward in secondary_indices - node_secondary = index_end - node_secondary_step = -1 - else - node_secondary = 1 - node_secondary_step = 1 - end - - for node in eachnode(dg) - # Get the normal direction on the primary element. - # Contravariant vectors at interfaces in negative coordinate direction - # are pointing inwards. This is handled by `get_normal_direction`. - normal_direction = get_normal_direction(primary_direction, - contravariant_vectors, - i_primary, j_primary, - primary_element) - - calc_interface_flux!(surface_flux_values, mesh, have_nonconservative_terms, - equations, - surface_integral, dg, cache, - interface, normal_direction, - node, primary_direction, primary_element, - node_secondary, secondary_direction, secondary_element) - - # Increment primary element indices to pull the normal direction - i_primary += i_primary_step - j_primary += j_primary_step - # Increment the surface node index along the secondary element - node_secondary += node_secondary_step - end - end - - return nothing -end -end # @muladd From 7e63a5a22645aee779c91f9bbfe1929274cd68ef Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Sun, 15 Mar 2026 12:58:02 +0100 Subject: [PATCH 20/46] comment --- src/solvers/dgsem_p4est/containers.jl | 2 +- src/solvers/dgsem_p4est/containers_parallel.jl | 3 +++ src/solvers/dgsem_t8code/containers.jl | 2 +- src/solvers/dgsem_t8code/containers_parallel.jl | 2 +- 4 files changed, 6 insertions(+), 3 deletions(-) diff --git a/src/solvers/dgsem_p4est/containers.jl b/src/solvers/dgsem_p4est/containers.jl index 9dea096b055..f109a4b4a9c 100644 --- a/src/solvers/dgsem_p4est/containers.jl +++ b/src/solvers/dgsem_p4est/containers.jl @@ -692,7 +692,7 @@ function reinitialize_containers!(mesh::P4estMesh, equations, dg::DGSEM, cache) # the number of iterations over the mesh in `p4est` init_surfaces!(interfaces, mortars, boundaries, mesh) - # Normal directions require that `node_indices` have been initialized + # init_normal_directions! requires that `node_indices` have been initialized init_normal_directions!(interfaces, dg.basis, elements) return nothing diff --git a/src/solvers/dgsem_p4est/containers_parallel.jl b/src/solvers/dgsem_p4est/containers_parallel.jl index e67ad26606d..f4cb9a31ede 100644 --- a/src/solvers/dgsem_p4est/containers_parallel.jl +++ b/src/solvers/dgsem_p4est/containers_parallel.jl @@ -314,6 +314,9 @@ function reinitialize_containers!(mesh::P4estMeshParallel, equations, dg::DGSEM, # the number of iterations over the mesh in p4est init_surfaces!(interfaces, mortars, boundaries, mpi_interfaces, mpi_mortars, mesh) + # init_normal_directions! requires that `node_indices` have been initialized + init_normal_directions!(interfaces, dg.basis, elements) + # re-initialize MPI cache @unpack mpi_cache = cache init_mpi_cache!(mpi_cache, mesh, mpi_interfaces, mpi_mortars, diff --git a/src/solvers/dgsem_t8code/containers.jl b/src/solvers/dgsem_t8code/containers.jl index ec5c6c9e9dd..1be31a1b20c 100644 --- a/src/solvers/dgsem_t8code/containers.jl +++ b/src/solvers/dgsem_t8code/containers.jl @@ -28,7 +28,7 @@ function reinitialize_containers!(mesh::T8codeMesh, equations, dg::DGSEM, cache) fill_mesh_info!(mesh, interfaces, mortars, boundaries, mesh.boundary_names) - # Normal directions require that `node_indices` have been initialized + # init_normal_directions! requires that `node_indices` have been initialized init_normal_directions!(interfaces, dg.basis, elements) return nothing diff --git a/src/solvers/dgsem_t8code/containers_parallel.jl b/src/solvers/dgsem_t8code/containers_parallel.jl index 27586c5b841..ac30a32daca 100644 --- a/src/solvers/dgsem_t8code/containers_parallel.jl +++ b/src/solvers/dgsem_t8code/containers_parallel.jl @@ -52,7 +52,7 @@ function reinitialize_containers!(mesh::T8codeMeshParallel, equations, dg::DGSEM init_normal_directions!(mpi_mortars, dg.basis, elements) exchange_normal_directions!(mpi_mortars, mpi_cache, mesh, nnodes(dg)) - # Normal directions require that `node_indices` have been initialized + # init_normal_directions! requires that `node_indices` have been initialized init_normal_directions!(interfaces, dg.basis, elements) return nothing From c9ff85fc7a74b76b5a2502213f916f1df7e3df16 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Sun, 15 Mar 2026 13:00:20 +0100 Subject: [PATCH 21/46] fmt --- src/solvers/dgsem_tree/dg_2d.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/solvers/dgsem_tree/dg_2d.jl b/src/solvers/dgsem_tree/dg_2d.jl index 69d1ffe3002..5e9b231151b 100644 --- a/src/solvers/dgsem_tree/dg_2d.jl +++ b/src/solvers/dgsem_tree/dg_2d.jl @@ -1271,7 +1271,8 @@ function calc_surface_integral!(du, u, end function calc_surface_integral!(du, u, - mesh::TreeMesh{2}, + mesh::Union{TreeMesh{2}, + StructuredMesh{2}, StructuredMeshView{2}}, equations, surface_integral::SurfaceIntegralWeakForm, dg::DGSEM{<:GaussLegendreBasis}, cache) @unpack boundary_interpolation_inverse_weights = dg.basis From 5f79f2e894d51e6b7c7a1ab341eb70aa41e07896 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Sun, 15 Mar 2026 13:53:11 +0100 Subject: [PATCH 22/46] other version t8coe --- src/solvers/dgsem_p4est/dg_2d.jl | 3 +- src/solvers/dgsem_t8code/containers.jl | 3 - .../dgsem_t8code/containers_parallel.jl | 3 - src/solvers/dgsem_t8code/dg_2d.jl | 72 +++++++++++++++++++ 4 files changed, 73 insertions(+), 8 deletions(-) create mode 100644 src/solvers/dgsem_t8code/dg_2d.jl diff --git a/src/solvers/dgsem_p4est/dg_2d.jl b/src/solvers/dgsem_p4est/dg_2d.jl index 544bcff4eef..b6be0e7a947 100644 --- a/src/solvers/dgsem_p4est/dg_2d.jl +++ b/src/solvers/dgsem_p4est/dg_2d.jl @@ -237,8 +237,7 @@ function prolong2interfaces!(cache, u, end function calc_interface_flux!(surface_flux_values, - mesh::Union{P4estMesh{2}, P4estMeshView{2}, - T8codeMesh{2}}, + mesh::Union{P4estMesh{2}, P4estMeshView{2}}, have_nonconservative_terms, equations, surface_integral, dg::DG, cache) @unpack neighbor_ids, node_indices, normal_directions = cache.interfaces diff --git a/src/solvers/dgsem_t8code/containers.jl b/src/solvers/dgsem_t8code/containers.jl index 1be31a1b20c..56db011f3b1 100644 --- a/src/solvers/dgsem_t8code/containers.jl +++ b/src/solvers/dgsem_t8code/containers.jl @@ -28,9 +28,6 @@ function reinitialize_containers!(mesh::T8codeMesh, equations, dg::DGSEM, cache) fill_mesh_info!(mesh, interfaces, mortars, boundaries, mesh.boundary_names) - # init_normal_directions! requires that `node_indices` have been initialized - init_normal_directions!(interfaces, dg.basis, elements) - return nothing end diff --git a/src/solvers/dgsem_t8code/containers_parallel.jl b/src/solvers/dgsem_t8code/containers_parallel.jl index ac30a32daca..65d4209a2c6 100644 --- a/src/solvers/dgsem_t8code/containers_parallel.jl +++ b/src/solvers/dgsem_t8code/containers_parallel.jl @@ -52,9 +52,6 @@ function reinitialize_containers!(mesh::T8codeMeshParallel, equations, dg::DGSEM init_normal_directions!(mpi_mortars, dg.basis, elements) exchange_normal_directions!(mpi_mortars, mpi_cache, mesh, nnodes(dg)) - # init_normal_directions! requires that `node_indices` have been initialized - init_normal_directions!(interfaces, dg.basis, elements) - return nothing end diff --git a/src/solvers/dgsem_t8code/dg_2d.jl b/src/solvers/dgsem_t8code/dg_2d.jl new file mode 100644 index 00000000000..308b268f196 --- /dev/null +++ b/src/solvers/dgsem_t8code/dg_2d.jl @@ -0,0 +1,72 @@ +@muladd begin +#! format: noindent + +# Use for `T8CodeMesh` the Gauss-Lobatto-Legendre (GLL) hard-coded variant, +# i.e., take interface normals from the outer volume nodes. +function calc_interface_flux!(surface_flux_values, + mesh::T8codeMesh{2}, + have_nonconservative_terms, + equations, surface_integral, dg::DG, cache) + @unpack neighbor_ids, node_indices = cache.interfaces + @unpack contravariant_vectors = cache.elements + index_range = eachnode(dg) + index_end = last(index_range) + + @threaded for interface in eachinterface(dg, cache) + # Get element and side index information on the primary element + primary_element = neighbor_ids[1, interface] + primary_indices = node_indices[1, interface] + primary_direction = indices2direction(primary_indices) + + # Create the local i,j indexing on the primary element used to pull normal direction information + i_primary_start, i_primary_step = index_to_start_step_2d(primary_indices[1], + index_range) + j_primary_start, j_primary_step = index_to_start_step_2d(primary_indices[2], + index_range) + + i_primary = i_primary_start + j_primary = j_primary_start + + # Get element and side index information on the secondary element + secondary_element = neighbor_ids[2, interface] + secondary_indices = node_indices[2, interface] + secondary_direction = indices2direction(secondary_indices) + + # Initiate the secondary index to be used in the surface for loop. + # This index on the primary side will always run forward but + # the secondary index might need to run backwards for flipped sides. + if :i_backward in secondary_indices + node_secondary = index_end + node_secondary_step = -1 + else + node_secondary = 1 + node_secondary_step = 1 + end + + for node in eachnode(dg) + # Get the normal direction on the primary element. + # Contravariant vectors at interfaces in negative coordinate direction + # are pointing inwards. This is handled by `get_normal_direction`. + normal_direction = get_normal_direction(primary_direction, + contravariant_vectors, + i_primary, j_primary, + primary_element) + + calc_interface_flux!(surface_flux_values, mesh, have_nonconservative_terms, + equations, + surface_integral, dg, cache, + interface, normal_direction, + node, primary_direction, primary_element, + node_secondary, secondary_direction, secondary_element) + + # Increment primary element indices to pull the normal direction + i_primary += i_primary_step + j_primary += j_primary_step + # Increment the surface node index along the secondary element + node_secondary += node_secondary_step + end + end + + return nothing +end +end # @muladd From 825872d7f585dcd61a5d4e9fd57c75bbbf80716f Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Sun, 15 Mar 2026 13:54:23 +0100 Subject: [PATCH 23/46] inc --- src/solvers/dgsem_t8code/dg.jl | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/solvers/dgsem_t8code/dg.jl b/src/solvers/dgsem_t8code/dg.jl index 47169c5d8e7..663044cdfac 100644 --- a/src/solvers/dgsem_t8code/dg.jl +++ b/src/solvers/dgsem_t8code/dg.jl @@ -34,5 +34,7 @@ include("containers_3d.jl") include("containers_parallel.jl") +include("dg_2d.jl") + include("dg_parallel.jl") end # @muladd From 022244d9f13cb92be917de936f50ee2a4654158e Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Sun, 15 Mar 2026 18:16:18 +0100 Subject: [PATCH 24/46] restore --- src/solvers/dgsem_p4est/dg_2d.jl | 73 ++++++++++++++++++++++++++++++- src/solvers/dgsem_t8code/dg.jl | 2 - src/solvers/dgsem_t8code/dg_2d.jl | 72 ------------------------------ 3 files changed, 72 insertions(+), 75 deletions(-) delete mode 100644 src/solvers/dgsem_t8code/dg_2d.jl diff --git a/src/solvers/dgsem_p4est/dg_2d.jl b/src/solvers/dgsem_p4est/dg_2d.jl index b6be0e7a947..449edab8fd6 100644 --- a/src/solvers/dgsem_p4est/dg_2d.jl +++ b/src/solvers/dgsem_p4est/dg_2d.jl @@ -236,10 +236,81 @@ function prolong2interfaces!(cache, u, return nothing end +# Take for Gauss-Lobatto-Legendre (GLL) the interface normals from the outer volume nodes. +function calc_interface_flux!(surface_flux_values, + mesh::Union{P4estMesh{2}, P4estMeshView{2}, + T8codeMesh{2}}, + have_nonconservative_terms, + equations, surface_integral, + dg::DGSEM{<:LobattoGaussLegendreBasis}, cache) + @unpack neighbor_ids, node_indices = cache.interfaces + @unpack contravariant_vectors = cache.elements + index_range = eachnode(dg) + index_end = last(index_range) + + @threaded for interface in eachinterface(dg, cache) + # Get element and side index information on the primary element + primary_element = neighbor_ids[1, interface] + primary_indices = node_indices[1, interface] + primary_direction = indices2direction(primary_indices) + + # Create the local i,j indexing on the primary element used to pull normal direction information + i_primary_start, i_primary_step = index_to_start_step_2d(primary_indices[1], + index_range) + j_primary_start, j_primary_step = index_to_start_step_2d(primary_indices[2], + index_range) + + i_primary = i_primary_start + j_primary = j_primary_start + + # Get element and side index information on the secondary element + secondary_element = neighbor_ids[2, interface] + secondary_indices = node_indices[2, interface] + secondary_direction = indices2direction(secondary_indices) + + # Initiate the secondary index to be used in the surface for loop. + # This index on the primary side will always run forward but + # the secondary index might need to run backwards for flipped sides. + if :i_backward in secondary_indices + node_secondary = index_end + node_secondary_step = -1 + else + node_secondary = 1 + node_secondary_step = 1 + end + + for node in eachnode(dg) + # Get the normal direction on the primary element. + # Contravariant vectors at interfaces in negative coordinate direction + # are pointing inwards. This is handled by `get_normal_direction`. + normal_direction = get_normal_direction(primary_direction, + contravariant_vectors, + i_primary, j_primary, + primary_element) + + calc_interface_flux!(surface_flux_values, mesh, have_nonconservative_terms, + equations, + surface_integral, dg, cache, + interface, normal_direction, + node, primary_direction, primary_element, + node_secondary, secondary_direction, secondary_element) + + # Increment primary element indices to pull the normal direction + i_primary += i_primary_step + j_primary += j_primary_step + # Increment the surface node index along the secondary element + node_secondary += node_secondary_step + end + end + + return nothing +end + function calc_interface_flux!(surface_flux_values, mesh::Union{P4estMesh{2}, P4estMeshView{2}}, have_nonconservative_terms, - equations, surface_integral, dg::DG, cache) + equations, surface_integral, + dg::DGSEM{<:GaussLegendreBasis}, cache) @unpack neighbor_ids, node_indices, normal_directions = cache.interfaces index_range = eachnode(dg) index_end = last(index_range) diff --git a/src/solvers/dgsem_t8code/dg.jl b/src/solvers/dgsem_t8code/dg.jl index 663044cdfac..47169c5d8e7 100644 --- a/src/solvers/dgsem_t8code/dg.jl +++ b/src/solvers/dgsem_t8code/dg.jl @@ -34,7 +34,5 @@ include("containers_3d.jl") include("containers_parallel.jl") -include("dg_2d.jl") - include("dg_parallel.jl") end # @muladd diff --git a/src/solvers/dgsem_t8code/dg_2d.jl b/src/solvers/dgsem_t8code/dg_2d.jl deleted file mode 100644 index 308b268f196..00000000000 --- a/src/solvers/dgsem_t8code/dg_2d.jl +++ /dev/null @@ -1,72 +0,0 @@ -@muladd begin -#! format: noindent - -# Use for `T8CodeMesh` the Gauss-Lobatto-Legendre (GLL) hard-coded variant, -# i.e., take interface normals from the outer volume nodes. -function calc_interface_flux!(surface_flux_values, - mesh::T8codeMesh{2}, - have_nonconservative_terms, - equations, surface_integral, dg::DG, cache) - @unpack neighbor_ids, node_indices = cache.interfaces - @unpack contravariant_vectors = cache.elements - index_range = eachnode(dg) - index_end = last(index_range) - - @threaded for interface in eachinterface(dg, cache) - # Get element and side index information on the primary element - primary_element = neighbor_ids[1, interface] - primary_indices = node_indices[1, interface] - primary_direction = indices2direction(primary_indices) - - # Create the local i,j indexing on the primary element used to pull normal direction information - i_primary_start, i_primary_step = index_to_start_step_2d(primary_indices[1], - index_range) - j_primary_start, j_primary_step = index_to_start_step_2d(primary_indices[2], - index_range) - - i_primary = i_primary_start - j_primary = j_primary_start - - # Get element and side index information on the secondary element - secondary_element = neighbor_ids[2, interface] - secondary_indices = node_indices[2, interface] - secondary_direction = indices2direction(secondary_indices) - - # Initiate the secondary index to be used in the surface for loop. - # This index on the primary side will always run forward but - # the secondary index might need to run backwards for flipped sides. - if :i_backward in secondary_indices - node_secondary = index_end - node_secondary_step = -1 - else - node_secondary = 1 - node_secondary_step = 1 - end - - for node in eachnode(dg) - # Get the normal direction on the primary element. - # Contravariant vectors at interfaces in negative coordinate direction - # are pointing inwards. This is handled by `get_normal_direction`. - normal_direction = get_normal_direction(primary_direction, - contravariant_vectors, - i_primary, j_primary, - primary_element) - - calc_interface_flux!(surface_flux_values, mesh, have_nonconservative_terms, - equations, - surface_integral, dg, cache, - interface, normal_direction, - node, primary_direction, primary_element, - node_secondary, secondary_direction, secondary_element) - - # Increment primary element indices to pull the normal direction - i_primary += i_primary_step - j_primary += j_primary_step - # Increment the surface node index along the secondary element - node_secondary += node_secondary_step - end - end - - return nothing -end -end # @muladd From 2c8291626d7d228528bee945bff8c78eecfd6a1a Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Sun, 15 Mar 2026 19:05:37 +0100 Subject: [PATCH 25/46] fix --- src/solvers/dgsem_p4est/dg_2d.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/solvers/dgsem_p4est/dg_2d.jl b/src/solvers/dgsem_p4est/dg_2d.jl index 449edab8fd6..6488eeb74f7 100644 --- a/src/solvers/dgsem_p4est/dg_2d.jl +++ b/src/solvers/dgsem_p4est/dg_2d.jl @@ -242,7 +242,7 @@ function calc_interface_flux!(surface_flux_values, T8codeMesh{2}}, have_nonconservative_terms, equations, surface_integral, - dg::DGSEM{<:LobattoGaussLegendreBasis}, cache) + dg::DGSEM{<:LobattoLegendreBasis}, cache) @unpack neighbor_ids, node_indices = cache.interfaces @unpack contravariant_vectors = cache.elements index_range = eachnode(dg) From bbf38c47cb1bf673f4bd5d6721ac1488c64120f1 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Tue, 17 Mar 2026 09:55:09 +0100 Subject: [PATCH 26/46] check if function is trixi atmo compatible --- src/solvers/dgsem_p4est/containers_2d.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/solvers/dgsem_p4est/containers_2d.jl b/src/solvers/dgsem_p4est/containers_2d.jl index 1934b44e84a..b2350d2305d 100644 --- a/src/solvers/dgsem_p4est/containers_2d.jl +++ b/src/solvers/dgsem_p4est/containers_2d.jl @@ -91,7 +91,7 @@ end # simply taken from the element face nodes. function init_normal_directions!(interfaces::P4estInterfaceContainer{2}, basis::LobattoLegendreBasis, - elements::P4estElementContainer{2}) + elements) @unpack neighbor_ids, node_indices, normal_directions = interfaces @unpack contravariant_vectors = elements index_range = eachnode(basis) @@ -131,7 +131,7 @@ end # computed from interpolation of the volume node normals to the surface nodes. function init_normal_directions!(interfaces::P4estInterfaceContainer{2}, basis::GaussLegendreBasis, - elements::P4estElementContainer{2}) + elements) @unpack neighbor_ids, node_indices, normal_directions = interfaces @unpack contravariant_vectors = elements @unpack boundary_interpolation = basis From ed1473eea01e5e1d1b38419395b7033042100f72 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Tue, 17 Mar 2026 10:14:42 +0100 Subject: [PATCH 27/46] restore --- src/solvers/dgsem_p4est/containers_2d.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/solvers/dgsem_p4est/containers_2d.jl b/src/solvers/dgsem_p4est/containers_2d.jl index b2350d2305d..1934b44e84a 100644 --- a/src/solvers/dgsem_p4est/containers_2d.jl +++ b/src/solvers/dgsem_p4est/containers_2d.jl @@ -91,7 +91,7 @@ end # simply taken from the element face nodes. function init_normal_directions!(interfaces::P4estInterfaceContainer{2}, basis::LobattoLegendreBasis, - elements) + elements::P4estElementContainer{2}) @unpack neighbor_ids, node_indices, normal_directions = interfaces @unpack contravariant_vectors = elements index_range = eachnode(basis) @@ -131,7 +131,7 @@ end # computed from interpolation of the volume node normals to the surface nodes. function init_normal_directions!(interfaces::P4estInterfaceContainer{2}, basis::GaussLegendreBasis, - elements) + elements::P4estElementContainer{2}) @unpack neighbor_ids, node_indices, normal_directions = interfaces @unpack contravariant_vectors = elements @unpack boundary_interpolation = basis From fb71bee6aef6649f5947abd097c50cf237c1b5b5 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Sun, 22 Mar 2026 12:37:58 +0100 Subject: [PATCH 28/46] comments --- src/meshes/p4est_mesh_view.jl | 5 +- src/solvers/dgsem_p4est/containers.jl | 64 +++++++++++++------ src/solvers/dgsem_p4est/containers_2d.jl | 40 ++---------- .../dgsem_p4est/containers_parallel_2d.jl | 2 + .../dgsem_p4est/containers_parallel_3d.jl | 2 + 5 files changed, 57 insertions(+), 56 deletions(-) diff --git a/src/meshes/p4est_mesh_view.jl b/src/meshes/p4est_mesh_view.jl index 5e2377dae7a..3d23a49f2bf 100644 --- a/src/meshes/p4est_mesh_view.jl +++ b/src/meshes/p4est_mesh_view.jl @@ -97,7 +97,10 @@ function extract_interfaces(mesh::P4estMeshView, interfaces_parent) # Copy relevant entries from parent mesh @views interfaces.u .= interfaces_parent.u[.., mask] - @views interfaces.normal_directions .= interfaces_parent.normal_directions[.., mask] + if interfaces.normal_directions !== nothing # Needed for Gauss-Legendre basis + @views interfaces.normal_directions .= interfaces_parent.normal_directions[.., + mask] + end @views interfaces.node_indices .= interfaces_parent.node_indices[.., mask] @views neighbor_ids = interfaces_parent.neighbor_ids[.., mask] diff --git a/src/solvers/dgsem_p4est/containers.jl b/src/solvers/dgsem_p4est/containers.jl index f109a4b4a9c..1e64844c88b 100644 --- a/src/solvers/dgsem_p4est/containers.jl +++ b/src/solvers/dgsem_p4est/containers.jl @@ -216,12 +216,12 @@ end mutable struct P4estInterfaceContainer{NDIMS, uEltype <: Real, RealT <: Real, NDIMSP1, NDIMSP2, uArray <: DenseArray{uEltype, NDIMSP2}, - NormalArray <: DenseArray{RealT, NDIMSP1}, + NormalArray, IdsMatrix <: DenseMatrix{Int}, IndicesMatrix <: DenseMatrix{NTuple{NDIMS, Symbol}}, uVector <: DenseVector{uEltype}, - NormalVector <: DenseVector{RealT}, + NormalVector, IdsVector <: DenseVector{Int}, IndicesVector <: DenseVector{NTuple{NDIMS, Symbol}}} <: @@ -238,6 +238,13 @@ mutable struct P4estInterfaceContainer{NDIMS, uEltype <: Real, RealT <: Real, _node_indices::IndicesVector end +# `trivial` means that the interface normals can be taken from +# the outer/surface element nodes +@inline trivial_interface_normals(::LobattoLegendreBasis) = true +# For Gauss-Legendre basis, the interface normals need to be interpolated, +# analogous to the surface values. +@inline trivial_interface_normals(::GaussLegendreBasis) = false + @inline function ninterfaces(interfaces::P4estInterfaceContainer) return size(interfaces.neighbor_ids, 2) end @@ -261,11 +268,17 @@ function Base.resize!(interfaces::P4estInterfaceContainer, capacity) (2, n_variables, ntuple(_ -> n_nodes, n_dims - 1)..., capacity)) - resize!(_normal_directions, n_dims * n_nodes^(n_dims - 1) * capacity) - interfaces.normal_directions = unsafe_wrap(ArrayType, pointer(_normal_directions), - (n_dims, - ntuple(_ -> n_nodes, n_dims - 1)..., - capacity)) + if _normal_directions === nothing # trivial interface normals, e.g. for Lobatto-Legendre basis + interfaces.normal_directions = nothing + else # non-trivial interface normals, e.g. for Gauss-Legendre basis + resize!(_normal_directions, n_dims * n_nodes^(n_dims - 1) * capacity) + interfaces.normal_directions = unsafe_wrap(ArrayType, + pointer(_normal_directions), + (n_dims, + ntuple(_ -> n_nodes, + n_dims - 1)..., + capacity)) + end resize!(_neighbor_ids, 2 * capacity) interfaces.neighbor_ids = unsafe_wrap(ArrayType, pointer(_neighbor_ids), @@ -295,13 +308,18 @@ function init_interfaces(mesh::Union{P4estMesh, P4estMeshView, T8codeMesh}, equa (2, nvariables(equations), ntuple(_ -> nnodes(basis), NDIMS - 1)..., n_interfaces)) - _normal_directions = Vector{RealT}(undef, - NDIMS * nnodes(basis)^(NDIMS - 1) * - n_interfaces) - normal_directions = unsafe_wrap(Array, pointer(_normal_directions), - (NDIMS, - ntuple(_ -> nnodes(basis), NDIMS - 1)..., - n_interfaces)) + if !trivial_interface_normals(basis) + _normal_directions = Vector{RealT}(undef, + NDIMS * nnodes(basis)^(NDIMS - 1) * + n_interfaces) + normal_directions = unsafe_wrap(Array, pointer(_normal_directions), + (NDIMS, + ntuple(_ -> nnodes(basis), NDIMS - 1)..., + n_interfaces)) + else + _normal_directions = nothing + normal_directions = nothing + end _neighbor_ids = Vector{Int}(undef, 2 * n_interfaces) neighbor_ids = unsafe_wrap(Array, pointer(_neighbor_ids), (2, n_interfaces)) @@ -343,25 +361,33 @@ function Adapt.parent_type(::Type{<:P4estInterfaceContainer{<:Any, <:Any, <:Any, end # Manual adapt_structure since we have aliasing memory -function Adapt.adapt_structure(to, interfaces::P4estInterfaceContainer) +function Adapt.adapt_structure(to, + interfaces::P4estInterfaceContainer{NDIMS, + uEltype, + RealT}) where { + NDIMS, + uEltype, + RealT + } # Adapt underlying storage _u = adapt(to, interfaces._u) - _normal_directions = adapt(to, interfaces._normal_directions) + _normal_directions = interfaces._normal_directions === nothing ? nothing : + adapt(to, interfaces._normal_directions) _neighbor_ids = adapt(to, interfaces._neighbor_ids) _node_indices = adapt(to, interfaces._node_indices) # Wrap arrays again u = unsafe_wrap_or_alloc(to, _u, size(interfaces.u)) - normal_directions = unsafe_wrap_or_alloc(to, _normal_directions, + normal_directions = _normal_directions === nothing ? nothing : # Check if interface normals are trivial + unsafe_wrap_or_alloc(to, _normal_directions, size(interfaces.normal_directions)) neighbor_ids = unsafe_wrap_or_alloc(to, _neighbor_ids, size(interfaces.neighbor_ids)) node_indices = unsafe_wrap_or_alloc(to, _node_indices, size(interfaces.node_indices)) - NDIMS = ndims(interfaces) new_type_params = (NDIMS, eltype(_u), - eltype(_normal_directions), + RealT, NDIMS + 1, NDIMS + 2, typeof(u), typeof(normal_directions), typeof(neighbor_ids), typeof(node_indices), diff --git a/src/solvers/dgsem_p4est/containers_2d.jl b/src/solvers/dgsem_p4est/containers_2d.jl index 1934b44e84a..e5ae8b5b4f5 100644 --- a/src/solvers/dgsem_p4est/containers_2d.jl +++ b/src/solvers/dgsem_p4est/containers_2d.jl @@ -87,44 +87,12 @@ function calc_node_coordinates!(node_coordinates, return node_coordinates end -# For Gauss-Lobatto-Legendre (GLL) nodes, the interface normals are -# simply taken from the element face nodes. +# For Gauss-Lobatto-Legendre (GLL) nodes, interface normals are computed on-the-fly +# during flux evaluation and do not need dedicated storage in the interface container. function init_normal_directions!(interfaces::P4estInterfaceContainer{2}, basis::LobattoLegendreBasis, elements::P4estElementContainer{2}) - @unpack neighbor_ids, node_indices, normal_directions = interfaces - @unpack contravariant_vectors = elements - index_range = eachnode(basis) - index_begin = first(index_range) - index_end = last(index_range) - - for interface in axes(neighbor_ids, 2) - primary_element = neighbor_ids[1, interface] - primary_indices = node_indices[1, interface] - primary_direction = indices2direction(primary_indices) - - i_primary, i_primary_step = index_to_start_step_2d(primary_indices[1], - index_begin, - index_end) - j_primary, j_primary_step = index_to_start_step_2d(primary_indices[2], - index_begin, - index_end) - - for i in index_range - # Retrieve normal direction at element/volume nodes - normal_direction = get_normal_direction(primary_direction, - contravariant_vectors, - i_primary, j_primary, - primary_element) - normal_directions[1, i, interface] = normal_direction[1] - normal_directions[2, i, interface] = normal_direction[2] - - i_primary += i_primary_step - j_primary += j_primary_step - end - end - - return interfaces + return nothing end # For Gauss-Legendre (GL) nodes, the interface normals are @@ -205,7 +173,7 @@ function init_normal_directions!(interfaces::P4estInterfaceContainer{2}, end end - return interfaces + return nothing end # Initialize node_indices of interface container diff --git a/src/solvers/dgsem_p4est/containers_parallel_2d.jl b/src/solvers/dgsem_p4est/containers_parallel_2d.jl index d531d33821b..31633cf2c88 100644 --- a/src/solvers/dgsem_p4est/containers_parallel_2d.jl +++ b/src/solvers/dgsem_p4est/containers_parallel_2d.jl @@ -77,5 +77,7 @@ function init_normal_directions!(mpi_mortars::P4estMPIMortarContainer{2}, basis, end end end + + return nothing end end # muladd diff --git a/src/solvers/dgsem_p4est/containers_parallel_3d.jl b/src/solvers/dgsem_p4est/containers_parallel_3d.jl index 56f0a543b97..b7419a2bebb 100644 --- a/src/solvers/dgsem_p4est/containers_parallel_3d.jl +++ b/src/solvers/dgsem_p4est/containers_parallel_3d.jl @@ -145,5 +145,7 @@ function init_normal_directions!(mpi_mortars::P4estMPIMortarContainer{3}, basis, end end end + + return nothing end end # muladd From 5e3c2e3345dbbb6711dbfc6e82ac4d053262b56a Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Sun, 22 Mar 2026 12:41:04 +0100 Subject: [PATCH 29/46] unions --- src/solvers/dgsem_p4est/containers.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/solvers/dgsem_p4est/containers.jl b/src/solvers/dgsem_p4est/containers.jl index 1e64844c88b..d03188b740a 100644 --- a/src/solvers/dgsem_p4est/containers.jl +++ b/src/solvers/dgsem_p4est/containers.jl @@ -216,12 +216,12 @@ end mutable struct P4estInterfaceContainer{NDIMS, uEltype <: Real, RealT <: Real, NDIMSP1, NDIMSP2, uArray <: DenseArray{uEltype, NDIMSP2}, - NormalArray, + NormalArray <: Union{DenseArray{RealT, NDIMSP1}, Nothing}, IdsMatrix <: DenseMatrix{Int}, IndicesMatrix <: DenseMatrix{NTuple{NDIMS, Symbol}}, uVector <: DenseVector{uEltype}, - NormalVector, + NormalVector <:Union{DenseVector{RealT}, Nothing}, IdsVector <: DenseVector{Int}, IndicesVector <: DenseVector{NTuple{NDIMS, Symbol}}} <: From 3f81bdfd083703a83739dc6b1f1ad842e694a3ac Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Sun, 22 Mar 2026 12:43:59 +0100 Subject: [PATCH 30/46] order yp eparams --- src/solvers/dgsem_p4est/containers.jl | 30 +++++++++++++-------------- 1 file changed, 15 insertions(+), 15 deletions(-) diff --git a/src/solvers/dgsem_p4est/containers.jl b/src/solvers/dgsem_p4est/containers.jl index d03188b740a..3621929f676 100644 --- a/src/solvers/dgsem_p4est/containers.jl +++ b/src/solvers/dgsem_p4est/containers.jl @@ -190,11 +190,8 @@ function Adapt.adapt_structure(to, size(elements.surface_flux_values)) new_type_params = (NDIMS, - RealT, - uEltype, - NDIMS + 1, - NDIMS + 2, - NDIMS + 3, + RealT, uEltype, + NDIMS + 1, NDIMS + 2, NDIMS + 3, typeof(inverse_jacobian), # ArrayRealTNDIMSP1 typeof(node_coordinates), # ArrayRealTNDIMSP2 typeof(jacobian_matrix), # ArrayRealTNDIMSP3 @@ -213,15 +210,17 @@ function Adapt.adapt_structure(to, _surface_flux_values) end -mutable struct P4estInterfaceContainer{NDIMS, uEltype <: Real, RealT <: Real, +mutable struct P4estInterfaceContainer{NDIMS, RealT <: Real, uEltype <: Real, NDIMSP1, NDIMSP2, uArray <: DenseArray{uEltype, NDIMSP2}, - NormalArray <: Union{DenseArray{RealT, NDIMSP1}, Nothing}, + NormalArray <: + Union{DenseArray{RealT, NDIMSP1}, Nothing}, IdsMatrix <: DenseMatrix{Int}, IndicesMatrix <: DenseMatrix{NTuple{NDIMS, Symbol}}, uVector <: DenseVector{uEltype}, - NormalVector <:Union{DenseVector{RealT}, Nothing}, + NormalVector <: + Union{DenseVector{RealT}, Nothing}, IdsVector <: DenseVector{Int}, IndicesVector <: DenseVector{NTuple{NDIMS, Symbol}}} <: @@ -249,8 +248,11 @@ end return size(interfaces.neighbor_ids, 2) end @inline Base.ndims(::P4estInterfaceContainer{NDIMS}) where {NDIMS} = NDIMS -@inline function Base.eltype(::P4estInterfaceContainer{NDIMS, uEltype}) where {NDIMS, - uEltype} +@inline function Base.eltype(::P4estInterfaceContainer{NDIMS, RealT, uEltype}) where { + NDIMS, + RealT, + uEltype + } return uEltype end @@ -327,7 +329,7 @@ function init_interfaces(mesh::Union{P4estMesh, P4estMeshView, T8codeMesh}, equa _node_indices = Vector{NTuple{NDIMS, Symbol}}(undef, 2 * n_interfaces) node_indices = unsafe_wrap(Array, pointer(_node_indices), (2, n_interfaces)) - interfaces = P4estInterfaceContainer{NDIMS, uEltype, RealT, + interfaces = P4estInterfaceContainer{NDIMS, RealT, uEltype, NDIMS + 1, NDIMS + 2, typeof(u), typeof(normal_directions), typeof(neighbor_ids), typeof(node_indices), @@ -386,8 +388,7 @@ function Adapt.adapt_structure(to, size(interfaces.node_indices)) new_type_params = (NDIMS, - eltype(_u), - RealT, + RealT, eltype(_u), NDIMS + 1, NDIMS + 2, typeof(u), typeof(normal_directions), typeof(neighbor_ids), typeof(node_indices), @@ -679,8 +680,7 @@ function Adapt.adapt_structure(to, mortars::P4estMortarContainer) NDIMS = ndims(mortars) new_type_params = (NDIMS, eltype(_u), - NDIMS + 1, - NDIMS + 3, + NDIMS + 1, NDIMS + 3, typeof(u), typeof(neighbor_ids), typeof(node_indices), typeof(_u), typeof(_neighbor_ids), typeof(_node_indices)) return P4estMortarContainer{new_type_params...}(u, neighbor_ids, node_indices, From 596854b374a61477bb19bddd46b097e75db2f919 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Sun, 22 Mar 2026 12:45:47 +0100 Subject: [PATCH 31/46] db --- src/solvers/dgsem_p4est/containers.jl | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/solvers/dgsem_p4est/containers.jl b/src/solvers/dgsem_p4est/containers.jl index 3621929f676..60162a868bc 100644 --- a/src/solvers/dgsem_p4est/containers.jl +++ b/src/solvers/dgsem_p4est/containers.jl @@ -365,12 +365,12 @@ end # Manual adapt_structure since we have aliasing memory function Adapt.adapt_structure(to, interfaces::P4estInterfaceContainer{NDIMS, - uEltype, - RealT}) where { - NDIMS, - uEltype, - RealT - } + RealT; + uEltype}) where { + NDIMS, + uEltype, + RealT + } # Adapt underlying storage _u = adapt(to, interfaces._u) _normal_directions = interfaces._normal_directions === nothing ? nothing : From 0e0104fb84aa61d9b99edbe03b64aa3a107d3297 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Sun, 22 Mar 2026 13:34:32 +0100 Subject: [PATCH 32/46] typo --- src/solvers/dgsem_p4est/containers.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/solvers/dgsem_p4est/containers.jl b/src/solvers/dgsem_p4est/containers.jl index 60162a868bc..bc7af295227 100644 --- a/src/solvers/dgsem_p4est/containers.jl +++ b/src/solvers/dgsem_p4est/containers.jl @@ -365,7 +365,7 @@ end # Manual adapt_structure since we have aliasing memory function Adapt.adapt_structure(to, interfaces::P4estInterfaceContainer{NDIMS, - RealT; + RealT, uEltype}) where { NDIMS, uEltype, From bd670cd5b1ca89c0c58ecc3f55e02b3238b29026 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Sun, 22 Mar 2026 20:36:09 +0100 Subject: [PATCH 33/46] lift --- src/solvers/dgsem_p4est/containers_2d.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/solvers/dgsem_p4est/containers_2d.jl b/src/solvers/dgsem_p4est/containers_2d.jl index e5ae8b5b4f5..3c46856f7cd 100644 --- a/src/solvers/dgsem_p4est/containers_2d.jl +++ b/src/solvers/dgsem_p4est/containers_2d.jl @@ -91,7 +91,7 @@ end # during flux evaluation and do not need dedicated storage in the interface container. function init_normal_directions!(interfaces::P4estInterfaceContainer{2}, basis::LobattoLegendreBasis, - elements::P4estElementContainer{2}) + elements) return nothing end @@ -99,7 +99,7 @@ end # computed from interpolation of the volume node normals to the surface nodes. function init_normal_directions!(interfaces::P4estInterfaceContainer{2}, basis::GaussLegendreBasis, - elements::P4estElementContainer{2}) + elements) @unpack neighbor_ids, node_indices, normal_directions = interfaces @unpack contravariant_vectors = elements @unpack boundary_interpolation = basis From 3263c11f987c0e6c9cfacd1482a276ddcaa21ded Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Sun, 22 Mar 2026 20:38:53 +0100 Subject: [PATCH 34/46] rev --- src/solvers/dgsem_p4est/containers_2d.jl | 6 ++---- src/solvers/dgsem_p4est/containers_3d.jl | 3 +-- src/solvers/dgsem_p4est/containers_parallel_2d.jl | 4 ++-- src/solvers/dgsem_p4est/containers_parallel_3d.jl | 4 ++-- 4 files changed, 7 insertions(+), 10 deletions(-) diff --git a/src/solvers/dgsem_p4est/containers_2d.jl b/src/solvers/dgsem_p4est/containers_2d.jl index 3c46856f7cd..0417f1b7212 100644 --- a/src/solvers/dgsem_p4est/containers_2d.jl +++ b/src/solvers/dgsem_p4est/containers_2d.jl @@ -90,16 +90,14 @@ end # For Gauss-Lobatto-Legendre (GLL) nodes, interface normals are computed on-the-fly # during flux evaluation and do not need dedicated storage in the interface container. function init_normal_directions!(interfaces::P4estInterfaceContainer{2}, - basis::LobattoLegendreBasis, - elements) + basis::LobattoLegendreBasis, elements) return nothing end # For Gauss-Legendre (GL) nodes, the interface normals are # computed from interpolation of the volume node normals to the surface nodes. function init_normal_directions!(interfaces::P4estInterfaceContainer{2}, - basis::GaussLegendreBasis, - elements) + basis::GaussLegendreBasis, elements) @unpack neighbor_ids, node_indices, normal_directions = interfaces @unpack contravariant_vectors = elements @unpack boundary_interpolation = basis diff --git a/src/solvers/dgsem_p4est/containers_3d.jl b/src/solvers/dgsem_p4est/containers_3d.jl index 55972861e64..c0e295153c7 100644 --- a/src/solvers/dgsem_p4est/containers_3d.jl +++ b/src/solvers/dgsem_p4est/containers_3d.jl @@ -78,8 +78,7 @@ end # Not yet implemented and needed for 3D function init_normal_directions!(interfaces::P4estInterfaceContainer{3}, - basis::LobattoLegendreBasis, - elements::P4estElementContainer{3}) + basis::LobattoLegendreBasis, elements) return nothing end diff --git a/src/solvers/dgsem_p4est/containers_parallel_2d.jl b/src/solvers/dgsem_p4est/containers_parallel_2d.jl index 31633cf2c88..003f97aab5f 100644 --- a/src/solvers/dgsem_p4est/containers_parallel_2d.jl +++ b/src/solvers/dgsem_p4est/containers_parallel_2d.jl @@ -39,8 +39,8 @@ end # Normal directions of small element surfaces are needed to calculate the mortar fluxes. Initialize # them for locally available small elements. -function init_normal_directions!(mpi_mortars::P4estMPIMortarContainer{2}, basis, - elements) +function init_normal_directions!(mpi_mortars::P4estMPIMortarContainer{2}, + basis, elements) @unpack local_neighbor_ids, local_neighbor_positions, node_indices = mpi_mortars @unpack contravariant_vectors = elements index_range = eachnode(basis) diff --git a/src/solvers/dgsem_p4est/containers_parallel_3d.jl b/src/solvers/dgsem_p4est/containers_parallel_3d.jl index b7419a2bebb..8a2c2939679 100644 --- a/src/solvers/dgsem_p4est/containers_parallel_3d.jl +++ b/src/solvers/dgsem_p4est/containers_parallel_3d.jl @@ -97,8 +97,8 @@ end # Normal directions of small element surfaces are needed to calculate the mortar fluxes. Initialize # them for locally available small elements. -function init_normal_directions!(mpi_mortars::P4estMPIMortarContainer{3}, basis, - elements) +function init_normal_directions!(mpi_mortars::P4estMPIMortarContainer{3}, + basis, elements) @unpack local_neighbor_ids, local_neighbor_positions, node_indices = mpi_mortars @unpack contravariant_vectors = elements index_range = eachnode(basis) From c30407718bfbb53c4070d6464c2d414db8b4c1e3 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Sun, 22 Mar 2026 20:38:57 +0100 Subject: [PATCH 35/46] ac --- examples/p4est_2d_dgsem/elixir_advection_flag_gauss_legendre.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/p4est_2d_dgsem/elixir_advection_flag_gauss_legendre.jl b/examples/p4est_2d_dgsem/elixir_advection_flag_gauss_legendre.jl index cbe37f087e3..ca2cf536684 100644 --- a/examples/p4est_2d_dgsem/elixir_advection_flag_gauss_legendre.jl +++ b/examples/p4est_2d_dgsem/elixir_advection_flag_gauss_legendre.jl @@ -38,7 +38,7 @@ summary_callback = SummaryCallback() analysis_callback = AnalysisCallback(semi, interval = 100) stepsize_callback = StepsizeCallback(cfl = 0.8) -callbacks = CallbackSet(summary_callback, #analysis_callback, +callbacks = CallbackSet(summary_callback, analysis_callback, stepsize_callback) ############################################################################### From 5fff783c55ad912760950c65972c7d25d8eb15f4 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Sun, 22 Mar 2026 23:11:28 +0100 Subject: [PATCH 36/46] rm --- src/solvers/dgsem_p4est/containers.jl | 7 ------- 1 file changed, 7 deletions(-) diff --git a/src/solvers/dgsem_p4est/containers.jl b/src/solvers/dgsem_p4est/containers.jl index bc7af295227..b92043428c1 100644 --- a/src/solvers/dgsem_p4est/containers.jl +++ b/src/solvers/dgsem_p4est/containers.jl @@ -248,13 +248,6 @@ end return size(interfaces.neighbor_ids, 2) end @inline Base.ndims(::P4estInterfaceContainer{NDIMS}) where {NDIMS} = NDIMS -@inline function Base.eltype(::P4estInterfaceContainer{NDIMS, RealT, uEltype}) where { - NDIMS, - RealT, - uEltype - } - return uEltype -end # See explanation of Base.resize! for the element container function Base.resize!(interfaces::P4estInterfaceContainer, capacity) From 9d598b13494efa1f085e793c8d7c44767e51e2c4 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Sat, 28 Mar 2026 14:16:11 +0100 Subject: [PATCH 37/46] adapt --- src/solvers/dgsem_p4est/dg_2d.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/solvers/dgsem_p4est/dg_2d.jl b/src/solvers/dgsem_p4est/dg_2d.jl index 6a4462323a7..f054c49498d 100644 --- a/src/solvers/dgsem_p4est/dg_2d.jl +++ b/src/solvers/dgsem_p4est/dg_2d.jl @@ -164,7 +164,7 @@ end return nothing end -function prolong2interfaces!(cache, u, +function prolong2interfaces!(backend::Nothing, cache, u, mesh::Union{P4estMesh{2}, P4estMeshView{2}}, equations, dg::DGSEM{<:GaussLegendreBasis}) @unpack interfaces = cache @@ -399,7 +399,7 @@ end return nothing end -function calc_interface_flux!(surface_flux_values, +function calc_interface_flux!(backend::Nothing, surface_flux_values, mesh::Union{P4estMesh{2}, P4estMeshView{2}}, have_nonconservative_terms, equations, surface_integral, @@ -1162,7 +1162,7 @@ end return nothing end -function calc_surface_integral!(du, u, +function calc_surface_integral!(backend::Nothing, du, u, mesh::Union{P4estMesh{2}, P4estMeshView{2}}, equations, surface_integral::SurfaceIntegralWeakForm, dg::DGSEM{<:GaussLegendreBasis}, cache) From 7047203be8868267b835bc556457fd49cad6278e Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Tue, 31 Mar 2026 23:24:11 +0200 Subject: [PATCH 38/46] adapt --- src/solvers/dgsem_p4est/dg_2d.jl | 111 +++++++++++++++++++------------ 1 file changed, 69 insertions(+), 42 deletions(-) diff --git a/src/solvers/dgsem_p4est/dg_2d.jl b/src/solvers/dgsem_p4est/dg_2d.jl index 002491d7ffc..b831129d832 100644 --- a/src/solvers/dgsem_p4est/dg_2d.jl +++ b/src/solvers/dgsem_p4est/dg_2d.jl @@ -272,7 +272,6 @@ function prolong2interfaces!(backend::Nothing, cache, u, return nothing end -# Take for Gauss-Lobatto-Legendre (GLL) the interface normals from the outer volume nodes. function calc_interface_flux!(backend::Nothing, surface_flux_values, mesh::Union{P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}}, @@ -280,6 +279,8 @@ function calc_interface_flux!(backend::Nothing, surface_flux_values, equations, surface_integral, dg::DGSEM{<:LobattoLegendreBasis}, cache) @unpack neighbor_ids, node_indices = cache.interfaces + # Take for Gauss-Lobatto-Legendre (GLL) the interface normals from the outer volume nodes, i.e., + # element data. @unpack contravariant_vectors = cache.elements index_range = eachnode(dg) @@ -299,7 +300,8 @@ function calc_interface_flux!(backend::Backend, surface_flux_values, mesh::Union{P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}}, have_nonconservative_terms, - equations, surface_integral, dg::DG, cache) + equations, surface_integral, + dg::DGSEM{<:LobattoLegendreBasis}, cache) ninterfaces(cache.interfaces) == 0 && return nothing @unpack neighbor_ids, node_indices = cache.interfaces @unpack contravariant_vectors = cache.elements @@ -337,7 +339,7 @@ end T8codeMesh{2}}}, have_nonconservative_terms, equations, surface_integral, - SolverT::Type{<:DG}, + SolverT::Type{<:DGSEM{<:LobattoLegendreBasis}}, u_interface, interface, neighbor_ids, node_indices, contravariant_vectors, @@ -385,8 +387,9 @@ end calc_interface_flux!(surface_flux_values, MeshT, have_nonconservative_terms, equations, surface_integral, SolverT, u_interface, - interface, normal_direction, node, primary_direction, - primary_element, node_secondary, + interface, normal_direction, node, + primary_direction, primary_element, + node_secondary, secondary_direction, secondary_element) # Increment primary element indices to pull the normal direction @@ -404,47 +407,71 @@ function calc_interface_flux!(backend::Nothing, surface_flux_values, have_nonconservative_terms, equations, surface_integral, dg::DGSEM{<:GaussLegendreBasis}, cache) - @unpack neighbor_ids, node_indices, normal_directions = cache.interfaces + @unpack neighbor_ids, node_indices = cache.interfaces + # Take for Gauss-Legendre (GL) the interface normals from the interfaces, i.e., + # interface data. + @unpack normal_directions = cache.interfaces index_range = eachnode(dg) - index_end = last(index_range) @threaded for interface in eachinterface(dg, cache) - # Get element and side index information on the primary element - primary_element = neighbor_ids[1, interface] - primary_indices = node_indices[1, interface] - primary_direction = indices2direction(primary_indices) - - # Get element and side index information on the secondary element - secondary_element = neighbor_ids[2, interface] - secondary_indices = node_indices[2, interface] - secondary_direction = indices2direction(secondary_indices) - - # Initiate the secondary index to be used in the surface for loop. - # This index on the primary side will always run forward but - # the secondary index might need to run backwards for flipped sides. - if :i_backward in secondary_indices - node_secondary = index_end - node_secondary_step = -1 - else - node_secondary = 1 - node_secondary_step = 1 - end + calc_interface_flux_per_interface!(surface_flux_values, typeof(mesh), + have_nonconservative_terms, + equations, surface_integral, typeof(dg), + cache.interfaces.u, interface, + neighbor_ids, node_indices, + normal_directions, index_range) + end - for node in eachnode(dg) - # This seems to be faster than `@views normal_direction = normal_directions[:, node, interface]` - normal_direction = SVector(normal_directions[1, node, interface], - normal_directions[2, node, interface]) - - calc_interface_flux!(surface_flux_values, mesh, have_nonconservative_terms, - equations, - surface_integral, dg, cache, - interface, normal_direction, - node, primary_direction, primary_element, - node_secondary, secondary_direction, secondary_element) - - # Increment the surface node index along the secondary element - node_secondary += node_secondary_step - end + return nothing +end + +@inline function calc_interface_flux_per_interface!(surface_flux_values, + MeshT::Type{<:Union{P4estMesh{2}, + P4estMeshView{2}}}, + have_nonconservative_terms, + equations, surface_integral, + SolverT::Type{<:DGSEM{<:GaussLegendreBasis}}, + u_interface, interface, + neighbor_ids, + node_indices, normal_directions, + index_range) + index_end = last(index_range) + + # Get element and side index information on the primary element + primary_element = neighbor_ids[1, interface] + primary_indices = node_indices[1, interface] + primary_direction = indices2direction(primary_indices) + + # Get element and side index information on the secondary element + secondary_element = neighbor_ids[2, interface] + secondary_indices = node_indices[2, interface] + secondary_direction = indices2direction(secondary_indices) + + # Initiate the secondary index to be used in the surface for loop. + # This index on the primary side will always run forward but + # the secondary index might need to run backwards for flipped sides. + if :i_backward in secondary_indices + node_secondary = index_end + node_secondary_step = -1 + else + node_secondary = 1 + node_secondary_step = 1 + end + + for node in index_range + # This seems to be faster than `@views normal_direction = normal_directions[:, node, interface]` + normal_direction = SVector(normal_directions[1, node, interface], + normal_directions[2, node, interface]) + + calc_interface_flux!(surface_flux_values, MeshT, have_nonconservative_terms, + equations, surface_integral, SolverT, u_interface, + interface, normal_direction, node, + primary_direction, primary_element, + node_secondary, + secondary_direction, secondary_element) + + # Increment the surface node index along the secondary element + node_secondary += node_secondary_step end return nothing From 96ff9fa65321da81e09cc15bfa97da5493e63f48 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Tue, 31 Mar 2026 23:31:36 +0200 Subject: [PATCH 39/46] cont --- src/solvers/dgsem_p4est/dg_2d.jl | 105 ++++++++++++++++++------------- 1 file changed, 60 insertions(+), 45 deletions(-) diff --git a/src/solvers/dgsem_p4est/dg_2d.jl b/src/solvers/dgsem_p4est/dg_2d.jl index b831129d832..acd6813547f 100644 --- a/src/solvers/dgsem_p4est/dg_2d.jl +++ b/src/solvers/dgsem_p4est/dg_2d.jl @@ -1104,7 +1104,7 @@ function calc_surface_integral!(backend::Nothing, du, u, mesh::Union{P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}}, equations, surface_integral::SurfaceIntegralWeakForm, - dg::DGSEM, cache) + dg::DGSEM{<:LobattoLegendreBasis}, cache) @unpack inverse_weights = dg.basis @unpack surface_flux_values = cache.elements @@ -1115,12 +1115,27 @@ function calc_surface_integral!(backend::Nothing, du, u, end end +function calc_surface_integral!(backend::Nothing, du, u, + mesh::Union{P4estMesh{2}, P4estMeshView{2}}, + equations, surface_integral::SurfaceIntegralWeakForm, + dg::DGSEM{<:GaussLegendreBasis}, cache) + @unpack boundary_interpolation_inverse_weights = dg.basis + @unpack surface_flux_values = cache.elements + + @threaded for element in eachelement(dg, cache) + calc_surface_integral_per_element!(du, typeof(mesh), equations, + surface_integral, dg, + boundary_interpolation_inverse_weights, + surface_flux_values, element) + end +end + function calc_surface_integral!(backend::Backend, du, u, mesh::Union{P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}}, equations, surface_integral::SurfaceIntegralWeakForm, - dg::DGSEM, cache) + dg::DGSEM{<:LobattoLegendreBasis}, cache) nelements(dg, cache) == 0 && return nothing @unpack inverse_weights = dg.basis @unpack surface_flux_values = cache.elements @@ -1137,7 +1152,8 @@ end T8codeMesh{2}}}, equations, surface_integral::SurfaceIntegralWeakForm, - dg::DGSEM, factor, + dg::DGSEM{<:LobattoLegendreBasis}, + factor, surface_flux_values) element = @index(Global) calc_surface_integral_per_element!(du, MeshT, equations, surface_integral, @@ -1150,9 +1166,9 @@ end T8codeMesh{2}}}, equations, surface_integral::SurfaceIntegralWeakForm, - dg::DGSEM, factor, - surface_flux_values, - element) + dg::DGSEM{<:LobattoLegendreBasis}, + factor, + surface_flux_values, element) # Note that all fluxes have been computed with outward-pointing normal vectors. # This computes the **negative** surface integral contribution, # i.e., M^{-1} * boundary_interpolation^T (which is for Gauss-Lobatto DGSEM just M^{-1} * B) @@ -1189,13 +1205,14 @@ end return nothing end -function calc_surface_integral!(backend::Nothing, du, u, - mesh::Union{P4estMesh{2}, P4estMeshView{2}}, - equations, surface_integral::SurfaceIntegralWeakForm, - dg::DGSEM{<:GaussLegendreBasis}, cache) - @unpack boundary_interpolation_inverse_weights = dg.basis - @unpack surface_flux_values = cache.elements - +function calc_surface_integral_per_element!(du, + ::Type{<:Union{P4estMesh{2}, + P4estMeshView{2}}}, + equations, + surface_integral::SurfaceIntegralWeakForm, + dg::DGSEM{<:GaussLegendreBasis}, + boundary_interpolation_inverse_weights, + surface_flux_values, element) # Note that all fluxes have been computed with outward-pointing normal vectors. # This computes the **negative** surface integral contribution, # i.e., M^{-1} * boundary_interpolation^T @@ -1203,39 +1220,37 @@ function calc_surface_integral!(backend::Nothing, du, u, # # We also use explicit assignments instead of `+=` to let `@muladd` turn these # into FMAs (see comment at the top of the file). - @threaded for element in eachelement(dg, cache) - for l in eachnode(dg) - for v in eachvariable(equations) - # Aliases for repeatedly accessed variables - surface_flux_minus_x = surface_flux_values[v, l, 1, element] - surface_flux_plus_x = surface_flux_values[v, l, 2, element] - for ii in eachnode(dg) - # surface at -x - du[v, ii, l, element] = (du[v, ii, l, element] + - surface_flux_minus_x * - boundary_interpolation_inverse_weights[ii, - 1]) - # surface at +x - du[v, ii, l, element] = (du[v, ii, l, element] + - surface_flux_plus_x * - boundary_interpolation_inverse_weights[ii, - 2]) - end + for l in eachnode(dg) + for v in eachvariable(equations) + # Aliases for repeatedly accessed variables + surface_flux_minus_x = surface_flux_values[v, l, 1, element] + surface_flux_plus_x = surface_flux_values[v, l, 2, element] + for ii in eachnode(dg) + # surface at -x + du[v, ii, l, element] = (du[v, ii, l, element] + + surface_flux_minus_x * + boundary_interpolation_inverse_weights[ii, + 1]) + # surface at +x + du[v, ii, l, element] = (du[v, ii, l, element] + + surface_flux_plus_x * + boundary_interpolation_inverse_weights[ii, + 2]) + end - surface_flux_minus_y = surface_flux_values[v, l, 3, element] - surface_flux_plus_y = surface_flux_values[v, l, 4, element] - for jj in eachnode(dg) - # surface at -y - du[v, l, jj, element] = (du[v, l, jj, element] + - surface_flux_minus_y * - boundary_interpolation_inverse_weights[jj, - 1]) - # surface at +y - du[v, l, jj, element] = (du[v, l, jj, element] + - surface_flux_plus_y * - boundary_interpolation_inverse_weights[jj, - 2]) - end + surface_flux_minus_y = surface_flux_values[v, l, 3, element] + surface_flux_plus_y = surface_flux_values[v, l, 4, element] + for jj in eachnode(dg) + # surface at -y + du[v, l, jj, element] = (du[v, l, jj, element] + + surface_flux_minus_y * + boundary_interpolation_inverse_weights[jj, + 1]) + # surface at +y + du[v, l, jj, element] = (du[v, l, jj, element] + + surface_flux_plus_y * + boundary_interpolation_inverse_weights[jj, + 2]) end end end From 69d898c48d77c9a8ddda23039b873eb70ef80335 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Wed, 1 Apr 2026 09:10:52 +0200 Subject: [PATCH 40/46] modularize --- src/solvers/dgsem_p4est/dg_2d.jl | 192 +++++++++++++++++-------------- 1 file changed, 104 insertions(+), 88 deletions(-) diff --git a/src/solvers/dgsem_p4est/dg_2d.jl b/src/solvers/dgsem_p4est/dg_2d.jl index acd6813547f..78051a01164 100644 --- a/src/solvers/dgsem_p4est/dg_2d.jl +++ b/src/solvers/dgsem_p4est/dg_2d.jl @@ -75,15 +75,15 @@ end function prolong2interfaces!(backend::Nothing, cache, u, mesh::Union{P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}}, - equations, dg::DG) + equations, dg::DGSEM{<:LobattoLegendreBasis}) @unpack interfaces = cache @unpack neighbor_ids, node_indices = cache.interfaces index_range = eachnode(dg) @threaded for interface in eachinterface(dg, cache) - prolong2interfaces_per_interface!(interfaces.u, u, interface, typeof(mesh), - equations, neighbor_ids, node_indices, - index_range) + prolong2interfaces_per_interface!(interfaces.u, u, interface, + typeof(mesh), equations, + neighbor_ids, node_indices, index_range) end return nothing end @@ -91,7 +91,7 @@ end function prolong2interfaces!(backend::Backend, cache, u, mesh::Union{P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}}, - equations, dg::DG) + equations, dg::DGSEM{<:LobattoLegendreBasis}) @unpack interfaces = cache ninterfaces(interfaces) == 0 && return nothing @unpack neighbor_ids, node_indices = cache.interfaces @@ -114,13 +114,13 @@ end neighbor_ids, node_indices, index_range) end +# Version for Gauss-Lobatto-Legendre @inline function prolong2interfaces_per_interface!(interfaces_u, u, interface, ::Type{<:Union{P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}}}, equations, neighbor_ids, - node_indices, - index_range) + node_indices, index_range) primary_element = neighbor_ids[1, interface] primary_indices = node_indices[1, interface] @@ -168,104 +168,120 @@ function prolong2interfaces!(backend::Nothing, cache, u, mesh::Union{P4estMesh{2}, P4estMeshView{2}}, equations, dg::DGSEM{<:GaussLegendreBasis}) @unpack interfaces = cache + @unpack neighbor_ids, node_indices = cache.interfaces @unpack boundary_interpolation = dg.basis index_range = eachnode(dg) @threaded for interface in eachinterface(dg, cache) - # Interpolate solution data from the primary element to the interface. - primary_element = interfaces.neighbor_ids[1, interface] - primary_indices = interfaces.node_indices[1, interface] + prolong2interfaces_per_interface!(interfaces.u, u, interface, + typeof(mesh), equations, + neighbor_ids, node_indices, index_range, + boundary_interpolation) + end - i_primary_start, i_primary_step = index_to_start_step_2d(primary_indices[1], - index_range) - j_primary_start, j_primary_step = index_to_start_step_2d(primary_indices[2], - index_range) - # The index direction is identified based on `{i,j}_{primary, secondary}_step`. - # For step = 0, the direction identified by this index is normal to the face. - # For step != 0 (1 or -1), the direction identified by this index is tangential to the face. - - # Note that in the current implementation, the interface will be - # "aligned at the primary element", i.e., the index of the primary side - # will always run forwards. - - i_primary = i_primary_start - j_primary = j_primary_start - - if i_primary_step == 0 - # i is the normal direction (constant), j varies along the surface - # => Interpolate in first/normal direction - # Interpolation side is governed by element orientation - side = interpolation_side(primary_indices[1]) - for i in eachnode(dg) - for v in eachvariable(equations) - u_primary = zero(eltype(interfaces.u)) - for ii in eachnode(dg) - u_primary = (u_primary + - u[v, ii, j_primary, primary_element] * - boundary_interpolation[ii, side]) - end - interfaces.u[1, v, i, interface] = u_primary + return nothing +end + +# Version for Gauss-Legendre, which requires passing in the boundary interpolation matrix +@inline function prolong2interfaces_per_interface!(interfaces_u, u, interface, + ::Type{<:Union{P4estMesh{2}, + P4estMeshView{2}}}, + equations, neighbor_ids, + node_indices, index_range, + boundary_interpolation) + # Interpolate solution data from the primary element to the interface. + primary_element = neighbor_ids[1, interface] + primary_indices = node_indices[1, interface] + + i_primary_start, i_primary_step = index_to_start_step_2d(primary_indices[1], + index_range) + j_primary_start, j_primary_step = index_to_start_step_2d(primary_indices[2], + index_range) + # The index direction is identified based on `{i,j}_{primary, secondary}_step`. + # For step = 0, the direction identified by this index is normal to the face. + # For step != 0 (1 or -1), the direction identified by this index is tangential to the face. + + # Note that in the current implementation, the interface will be + # "aligned at the primary element", i.e., the index of the primary side + # will always run forwards. + + i_primary = i_primary_start + j_primary = j_primary_start + + if i_primary_step == 0 + # i is the normal direction (constant), j varies along the surface + # => Interpolate in first/normal direction + # Interpolation side is governed by element orientation + side = interpolation_side(primary_indices[1]) + for i in index_range + for v in eachvariable(equations) + u_primary = zero(eltype(interfaces_u)) + for ii in index_range + u_primary = (u_primary + + u[v, ii, j_primary, primary_element] * + boundary_interpolation[ii, side]) end - j_primary += j_primary_step # incrementing j_primary suffices + interfaces_u[1, v, i, interface] = u_primary end - else # j_primary_step == 0 - # j is the normal direction (constant), i varies along the surface - # => Interpolate in second/normal direction - # Interpolation side is governed by element orientation - side = interpolation_side(primary_indices[2]) - for i in eachnode(dg) - for v in eachvariable(equations) - u_primary = zero(eltype(interfaces.u)) - for jj in eachnode(dg) - u_primary = (u_primary + - u[v, i_primary, jj, primary_element] * - boundary_interpolation[jj, side]) - end - interfaces.u[1, v, i, interface] = u_primary + j_primary += j_primary_step # incrementing j_primary suffices + end + else # j_primary_step == 0 + # j is the normal direction (constant), i varies along the surface + # => Interpolate in second/normal direction + # Interpolation side is governed by element orientation + side = interpolation_side(primary_indices[2]) + for i in index_range + for v in eachvariable(equations) + u_primary = zero(eltype(interfaces_u)) + for jj in index_range + u_primary = (u_primary + + u[v, i_primary, jj, primary_element] * + boundary_interpolation[jj, side]) end - i_primary += i_primary_step # incrementing i_primary suffices + interfaces_u[1, v, i, interface] = u_primary end + i_primary += i_primary_step # incrementing i_primary suffices end + end - # Interpolate solution data from the secondary element to the interface. - secondary_element = interfaces.neighbor_ids[2, interface] - secondary_indices = interfaces.node_indices[2, interface] + # Interpolate solution data from the secondary element to the interface. + secondary_element = neighbor_ids[2, interface] + secondary_indices = node_indices[2, interface] - i_secondary_start, i_secondary_step = index_to_start_step_2d(secondary_indices[1], - index_range) - j_secondary_start, j_secondary_step = index_to_start_step_2d(secondary_indices[2], - index_range) + i_secondary_start, i_secondary_step = index_to_start_step_2d(secondary_indices[1], + index_range) + j_secondary_start, j_secondary_step = index_to_start_step_2d(secondary_indices[2], + index_range) - i_secondary = i_secondary_start - j_secondary = j_secondary_start - if i_secondary_step == 0 - side = interpolation_side(secondary_indices[1]) - for i in eachnode(dg) - for v in eachvariable(equations) - u_secondary = zero(eltype(interfaces.u)) - for ii in eachnode(dg) - u_secondary = (u_secondary + - u[v, ii, j_secondary, secondary_element] * - boundary_interpolation[ii, side]) - end - interfaces.u[2, v, i, interface] = u_secondary + i_secondary = i_secondary_start + j_secondary = j_secondary_start + if i_secondary_step == 0 + side = interpolation_side(secondary_indices[1]) + for i in index_range + for v in eachvariable(equations) + u_secondary = zero(eltype(interfaces_u)) + for ii in index_range + u_secondary = (u_secondary + + u[v, ii, j_secondary, secondary_element] * + boundary_interpolation[ii, side]) end - j_secondary += j_secondary_step # incrementing j_secondary suffices + interfaces_u[2, v, i, interface] = u_secondary end - else # j_secondary_step == 0 - side = interpolation_side(secondary_indices[2]) - for i in eachnode(dg) - for v in eachvariable(equations) - u_secondary = zero(eltype(interfaces.u)) - for jj in eachnode(dg) - u_secondary = (u_secondary + - u[v, i_secondary, jj, secondary_element] * - boundary_interpolation[jj, side]) - end - interfaces.u[2, v, i, interface] = u_secondary + j_secondary += j_secondary_step # incrementing j_secondary suffices + end + else # j_secondary_step == 0 + side = interpolation_side(secondary_indices[2]) + for i in index_range + for v in eachvariable(equations) + u_secondary = zero(eltype(interfaces_u)) + for jj in index_range + u_secondary = (u_secondary + + u[v, i_secondary, jj, secondary_element] * + boundary_interpolation[jj, side]) end - i_secondary += i_secondary_step # incrementing i_secondary suffices + interfaces_u[2, v, i, interface] = u_secondary end + i_secondary += i_secondary_step # incrementing i_secondary suffices end end From 6c80410d69cc9e6c41a161a68d7a44e094e75cf9 Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Thu, 2 Apr 2026 12:03:48 +0200 Subject: [PATCH 41/46] Update examples/p4est_2d_dgsem/elixir_advection_flag_gauss_legendre.jl Co-authored-by: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com> --- .../p4est_2d_dgsem/elixir_advection_flag_gauss_legendre.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/p4est_2d_dgsem/elixir_advection_flag_gauss_legendre.jl b/examples/p4est_2d_dgsem/elixir_advection_flag_gauss_legendre.jl index ca2cf536684..0b09bdc5606 100644 --- a/examples/p4est_2d_dgsem/elixir_advection_flag_gauss_legendre.jl +++ b/examples/p4est_2d_dgsem/elixir_advection_flag_gauss_legendre.jl @@ -14,8 +14,8 @@ solver = DGSEM(polydeg = 4, surface_flux = flux_lax_friedrichs, # lower and upper faces are sinus curves, left and right are vertical lines. f1(s) = SVector(-1.0, s - 1.0) f2(s) = SVector(1.0, s + 1.0) -f3(s) = SVector(s, -1.0 + sin(0.5 * pi * s)) -f4(s) = SVector(s, 1.0 + sin(0.5 * pi * s)) +f3(s) = SVector(s, -1.0 + sinpi(0.5 * s)) +f4(s) = SVector(s, 1.0 + sinpi(0.5 * s)) # Create P4estMesh with 3 x 2 trees and 6 x 4 elements, # approximate the geometry with a smaller polydeg for testing. From 4fe7262efff3d70fadb220871a5ee56a3fbfb609 Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Mon, 6 Apr 2026 09:16:08 +0200 Subject: [PATCH 42/46] Update src/solvers/dgsem_p4est/containers_2d.jl Co-authored-by: Hendrik Ranocha --- src/solvers/dgsem_p4est/containers_2d.jl | 8 ++------ 1 file changed, 2 insertions(+), 6 deletions(-) diff --git a/src/solvers/dgsem_p4est/containers_2d.jl b/src/solvers/dgsem_p4est/containers_2d.jl index 0417f1b7212..a0f7accf6c3 100644 --- a/src/solvers/dgsem_p4est/containers_2d.jl +++ b/src/solvers/dgsem_p4est/containers_2d.jl @@ -102,8 +102,6 @@ function init_normal_directions!(interfaces::P4estInterfaceContainer{2}, @unpack contravariant_vectors = elements @unpack boundary_interpolation = basis index_range = eachnode(basis) - index_begin = first(index_range) - index_end = last(index_range) for interface in axes(neighbor_ids, 2) primary_element = neighbor_ids[1, interface] @@ -111,11 +109,9 @@ function init_normal_directions!(interfaces::P4estInterfaceContainer{2}, primary_direction = indices2direction(primary_indices) i_primary_start, i_primary_step = index_to_start_step_2d(primary_indices[1], - index_begin, - index_end) + index_range) j_primary_start, j_primary_step = index_to_start_step_2d(primary_indices[2], - index_begin, - index_end) + index_range) # The index direction is identified based on `{i,j}_{primary, secondary}_step`. # For step = 0, the direction identified by this index is normal to the face. From e2c00299b39dc5eee8218c179a6f7254aa8fea7d Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Mon, 6 Apr 2026 09:17:29 +0200 Subject: [PATCH 43/46] rev --- src/solvers/dgsem_p4est/dg_2d.jl | 11 ++++------- 1 file changed, 4 insertions(+), 7 deletions(-) diff --git a/src/solvers/dgsem_p4est/dg_2d.jl b/src/solvers/dgsem_p4est/dg_2d.jl index 78051a01164..9af6029aca4 100644 --- a/src/solvers/dgsem_p4est/dg_2d.jl +++ b/src/solvers/dgsem_p4est/dg_2d.jl @@ -48,7 +48,10 @@ end # i_volume += i_volume_step # j_volume += j_volume_step # end -@inline function index_to_start_step_2d(index::Symbol, index_begin, index_end) +@inline function index_to_start_step_2d(index::Symbol, index_range) + index_begin = first(index_range) + index_end = last(index_range) + if index === :begin return index_begin, 0 elseif index === :end @@ -59,12 +62,6 @@ end return index_end, -1 end end -@inline function index_to_start_step_2d(index::Symbol, index_range) - index_begin = first(index_range) - index_end = last(index_range) - - return index_to_start_step_2d(index, index_begin, index_end) -end # Infer interpolation side, i.e., left (1) or right (2) for an element. # Required for boundary interpolation with Gauss-Legendre nodes. From ad3554230e2e14ab60cddf26b1995bd0d023528f Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Mon, 6 Apr 2026 09:22:06 +0200 Subject: [PATCH 44/46] rm --- src/solvers/dgsem_p4est/dg_2d.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/src/solvers/dgsem_p4est/dg_2d.jl b/src/solvers/dgsem_p4est/dg_2d.jl index 9af6029aca4..f60ece112c1 100644 --- a/src/solvers/dgsem_p4est/dg_2d.jl +++ b/src/solvers/dgsem_p4est/dg_2d.jl @@ -26,7 +26,6 @@ function create_cache(mesh::Union{P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}} return cache end -# index_to_start_step_2d(index::Symbol, index_begin, index_end) # index_to_start_step_2d(index::Symbol, index_range) # # Given a symbolic `index` and an `indexrange` (usually `eachnode(dg)`), From 46080710dea05d769481b604cc6f6208d2ce438d Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Mon, 6 Apr 2026 09:25:45 +0200 Subject: [PATCH 45/46] Update examples/p4est_2d_dgsem/elixir_advection_flag_gauss_legendre.jl Co-authored-by: Hendrik Ranocha --- examples/p4est_2d_dgsem/elixir_advection_flag_gauss_legendre.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/p4est_2d_dgsem/elixir_advection_flag_gauss_legendre.jl b/examples/p4est_2d_dgsem/elixir_advection_flag_gauss_legendre.jl index 0b09bdc5606..815317ae26c 100644 --- a/examples/p4est_2d_dgsem/elixir_advection_flag_gauss_legendre.jl +++ b/examples/p4est_2d_dgsem/elixir_advection_flag_gauss_legendre.jl @@ -36,7 +36,7 @@ ode = semidiscretize(semi, (0.0, 0.2)) summary_callback = SummaryCallback() analysis_callback = AnalysisCallback(semi, interval = 100) -stepsize_callback = StepsizeCallback(cfl = 0.8) +stepsize_callback = StepsizeCallback(cfl = 0.5) callbacks = CallbackSet(summary_callback, analysis_callback, stepsize_callback) From 40acf2bb5e0570c656bc0b9714ad435d6a2f514a Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Mon, 6 Apr 2026 09:36:24 +0200 Subject: [PATCH 46/46] test --- test/test_p4est_2d.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_p4est_2d.jl b/test/test_p4est_2d.jl index 41e6a13bb61..5147ce29aa3 100644 --- a/test/test_p4est_2d.jl +++ b/test/test_p4est_2d.jl @@ -76,7 +76,7 @@ end @trixi_testset "elixir_advection_flag_gauss_legendre.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_flag_gauss_legendre.jl"), - l2=[0.00047332539212915036], linf=[0.0021788728950808967]) + l2=[0.0004734270965231037], linf=[0.002206239481024719]) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) @test_allocations(Trixi.rhs!, semi, sol, 1000)