Refactor StructuredMesh

src/semidiscretization/semidiscretization_coupled.jl

@@ -646,7 +646,7 @@ end
 ### DGSEM/structured
 ################################################################################
 
-@inline function calc_boundary_flux_by_direction!(surface_flux_values, u, t,
+@inline function calc_boundary_flux_by_direction!(surface_flux_values, t,
                                                   orientation,
                                                   boundary_condition::BoundaryConditionCoupled,
                                                   mesh::Union{StructuredMesh,
@@ -656,12 +656,15 @@ end
                                                   surface_integral, dg::DG, cache,
                                                   direction, node_indices,
                                                   surface_node_indices, element)
-    @unpack node_coordinates, contravariant_vectors, inverse_jacobian = cache.elements
+    @unpack node_coordinates, contravariant_vectors, inverse_jacobian, interfaces_u = cache.elements
+    # Boundary values are for `StructuredMesh` stored in the interface datastructure
+    boundaries_u = interfaces_u
     @unpack surface_flux = surface_integral
 
     cell_indices = get_boundary_indices(element, orientation, mesh)
 
-    u_inner = get_node_vars(u, equations, dg, node_indices..., element)
+    u_inner = get_node_vars(boundaries_u, equations, dg, surface_node_indices...,
+                            direction, element)
 
     # If the mapping is orientation-reversing, the contravariant vectors' orientation
     # is reversed as well. The normal vector must be oriented in the direction
@@ -686,7 +689,7 @@ end
     return nothing
 end
 
-@inline function calc_boundary_flux_by_direction!(surface_flux_values, u, t,
+@inline function calc_boundary_flux_by_direction!(surface_flux_values, t,
                                                   orientation,
                                                   boundary_condition::BoundaryConditionCoupled,
                                                   mesh::Union{StructuredMesh,
@@ -696,12 +699,15 @@ end
                                                   surface_integral, dg::DG, cache,
                                                   direction, node_indices,
                                                   surface_node_indices, element)
-    @unpack node_coordinates, contravariant_vectors, inverse_jacobian = cache.elements
+    @unpack node_coordinates, contravariant_vectors, inverse_jacobian, interfaces_u = cache.elements
+    # Boundary values are for `StructuredMesh` stored in the interface datastructure
+    boundaries_u = interfaces_u
     @unpack surface_flux = surface_integral
 
     cell_indices = get_boundary_indices(element, orientation, mesh)
 
-    u_inner = get_node_vars(u, equations, dg, node_indices..., element)
+    u_inner = get_node_vars(boundaries_u, equations, dg, surface_node_indices...,
+                            direction, element)
 
     # If the mapping is orientation-reversing, the contravariant vectors' orientation
     # is reversed as well. The normal vector must be oriented in the direction

src/solvers/dgsem_structured/containers.jl

@@ -22,6 +22,9 @@ struct StructuredElementContainer{NDIMS, RealT <: Real, uEltype <: Real,
     # 1/J where J is the Jacobian determinant (determinant of Jacobian matrix)
     inverse_jacobian::Array{RealT, NDIMSP1} # [node_i, node_j, node_k, element]
 
+    # Buffer for solution values at interfaces (filled by `prolong2interfaces!`)
+    interfaces_u::Array{uEltype, NDIMSP2} # [variable, i, j, direction, element]
+
     # Buffer for calculated surface flux
     surface_flux_values::Array{uEltype, NDIMSP2} # [variable, i, j, direction, element]
 end
@@ -44,6 +47,10 @@ function init_elements(mesh::Union{StructuredMesh{NDIMS, RealT},
     inverse_jacobian = Array{RealT, NDIMS + 1}(undef,
                                                ntuple(_ -> nnodes(basis), NDIMS)...,
                                                nelements)
+    interfaces_u = Array{uEltype, NDIMS + 2}(undef, nvariables(equations),
+                                             ntuple(_ -> nnodes(basis),
+                                                    NDIMS - 1)..., NDIMS * 2,
+                                             nelements)
     surface_flux_values = Array{uEltype, NDIMS + 2}(undef, nvariables(equations),
                                                     ntuple(_ -> nnodes(basis),
                                                            NDIMS - 1)..., NDIMS * 2,
@@ -55,6 +62,7 @@ function init_elements(mesh::Union{StructuredMesh{NDIMS, RealT},
                                                                            jacobian_matrix,
                                                                            contravariant_vectors,
                                                                            inverse_jacobian,
+                                                                           interfaces_u,
                                                                            surface_flux_values)
 
     init_elements!(elements, mesh, basis)

src/solvers/dgsem_structured/dg.jl

@@ -30,7 +30,7 @@ end
 end
 
 # Dimension agnostic, i.e., valid for all 1D, 2D, and 3D `StructuredMesh`es.
-function calc_boundary_flux!(cache, u, t, boundary_condition::BoundaryConditionPeriodic,
+function calc_boundary_flux!(cache, t, boundary_condition::BoundaryConditionPeriodic,
                              mesh::StructuredMesh, equations, surface_integral,
                              dg::DG)
     @assert isperiodic(mesh)
@@ -52,16 +52,25 @@ function rhs!(du, u, t,
                               dg.volume_integral, dg, cache)
     end
 
-    # Calculate interface and boundary fluxes
+    # Prolong solution to interfaces
+    @trixi_timeit timer() "prolong2interfaces" begin
+        prolong2interfaces!(cache, u, mesh, equations, dg)
+    end
+
+    # Calculate interface fluxes
     @trixi_timeit timer() "interface flux" begin
-        calc_interface_flux!(cache, u, mesh,
+        calc_interface_flux!(cache.elements.surface_flux_values, mesh,
                              have_nonconservative_terms(equations), equations,
-                             dg.surface_integral, dg)
+                             dg.surface_integral, dg, cache)
     end
 
+    # `prolong2boundaries!` is not required for `StructuredMesh` since boundary values 
+    # are stored in the interface datastructure (`interfaces_u`),
+    # so we can directly calculate the boundary fluxes without prolongation.
+
     # Calculate boundary fluxes
     @trixi_timeit timer() "boundary flux" begin
-        calc_boundary_flux!(cache, u, t, boundary_conditions, mesh, equations,
+        calc_boundary_flux!(cache, t, boundary_conditions, mesh, equations,
                             dg.surface_integral, dg)
     end
 
@@ -82,7 +91,7 @@ function rhs!(du, u, t,
     return nothing
 end
 
-@inline function calc_boundary_flux_by_direction!(surface_flux_values, u, t,
+@inline function calc_boundary_flux_by_direction!(surface_flux_values, t,
                                                   orientation,
                                                   boundary_condition::BoundaryConditionPeriodic,
                                                   mesh::Union{StructuredMesh,
@@ -96,7 +105,7 @@ end
     return nothing
 end
 
-@inline function calc_boundary_flux_by_direction!(surface_flux_values, u, t,
+@inline function calc_boundary_flux_by_direction!(surface_flux_values, t,
                                                   orientation,
                                                   boundary_condition::BoundaryConditionPeriodic,
                                                   mesh::Union{StructuredMesh,
@@ -110,7 +119,7 @@ end
     return nothing
 end
 
-@inline function calc_boundary_flux_by_direction!(surface_flux_values, u, t,
+@inline function calc_boundary_flux_by_direction!(surface_flux_values, t,
                                                   orientation,
                                                   boundary_condition,
                                                   mesh::Union{StructuredMesh,
@@ -120,10 +129,13 @@ end
                                                   surface_integral, dg::DG, cache,
                                                   direction, node_indices,
                                                   surface_node_indices, element)
-    @unpack node_coordinates, contravariant_vectors, inverse_jacobian = cache.elements
+    @unpack node_coordinates, contravariant_vectors, inverse_jacobian, interfaces_u = cache.elements
+    # Boundary values are for `StructuredMesh` stored in the interface datastructure
+    boundaries_u = interfaces_u
     @unpack surface_flux = surface_integral
 
-    u_inner = get_node_vars(u, equations, dg, node_indices..., element)
+    u_inner = get_node_vars(boundaries_u, equations, dg, surface_node_indices...,
+                            direction, element)
     x = get_node_coords(node_coordinates, equations, dg, node_indices..., element)
 
     # If the mapping is orientation-reversing, the contravariant vectors' orientation
@@ -142,14 +154,15 @@ end
     flux = sign_jacobian *
            boundary_condition(u_inner, normal, direction, x, t, surface_flux, equations)
 
+    # Only flux contribution for boundary element, boundary face is the boundary flux
     for v in eachvariable(equations)
         surface_flux_values[v, surface_node_indices..., direction, element] = flux[v]
     end
 
     return nothing
 end
 
-@inline function calc_boundary_flux_by_direction!(surface_flux_values, u, t,
+@inline function calc_boundary_flux_by_direction!(surface_flux_values, t,
                                                   orientation,
                                                   boundary_condition,
                                                   mesh::Union{StructuredMesh,
@@ -159,10 +172,13 @@ end
                                                   surface_integral, dg::DG, cache,
                                                   direction, node_indices,
                                                   surface_node_indices, element)
-    @unpack node_coordinates, contravariant_vectors, inverse_jacobian = cache.elements
+    @unpack node_coordinates, contravariant_vectors, inverse_jacobian, interfaces_u = cache.elements
+    # Boundary values are for `StructuredMesh` stored in the interface datastructure
+    boundaries_u = interfaces_u
     @unpack surface_flux = surface_integral
 
-    u_inner = get_node_vars(u, equations, dg, node_indices..., element)
+    u_inner = get_node_vars(boundaries_u, equations, dg, surface_node_indices...,
+                            direction, element)
     x = get_node_coords(node_coordinates, equations, dg, node_indices..., element)
 
     # If the mapping is orientation-reversing, the contravariant vectors' orientation
@@ -181,6 +197,7 @@ end
     flux, noncons_flux = boundary_condition(u_inner, normal, direction, x, t,
                                             surface_flux, equations)
 
+    # Only flux contribution for boundary element, boundary face is the boundary flux
     for v in eachvariable(equations)
         surface_flux_values[v, surface_node_indices..., direction, element] = sign_jacobian *
                                                                               (flux[v] +

src/solvers/dgsem_structured/dg_1d.jl

@@ -5,63 +5,84 @@
 @muladd begin
 #! format: noindent
 
-function calc_interface_flux!(cache, u, mesh::StructuredMesh{1},
+function prolong2interfaces!(cache, u, mesh::StructuredMesh{1}, equations, dg::DG)
+    @unpack interfaces_u = cache.elements
+
+    @threaded for element in eachelement(dg, cache)
+        # Negative side (direction 1, left/negative x face)
+        for v in eachvariable(equations)
+            interfaces_u[v, 1, element] = u[v, 1, element]
+        end
+        # Positive side (direction 2, right/positive x face)
+        for v in eachvariable(equations)
+            interfaces_u[v, 2, element] = u[v, nnodes(dg), element]
+        end
+    end
+
+    return nothing
+end
+
+function calc_interface_flux!(surface_flux_values, mesh::StructuredMesh{1},
                               nonconservative_terms, # can be True/False
-                              equations, surface_integral, dg::DG)
+                              equations, surface_integral, dg::DG, cache)
     @unpack surface_flux = surface_integral
+    @unpack interfaces_u = cache.elements
 
     @threaded for element in eachelement(dg, cache)
         left_element = cache.elements.left_neighbors[1, element]
         # => `element` is the right element of the interface
 
         if left_element > 0 # left_element = 0 at boundaries
-            u_ll = get_node_vars(u, equations, dg, nnodes(dg), left_element)
-            u_rr = get_node_vars(u, equations, dg, 1, element)
+            u_ll = get_node_vars(interfaces_u, equations, dg, 2, left_element)
+            u_rr = get_node_vars(interfaces_u, equations, dg, 1, element)
 
             f1 = surface_flux(u_ll, u_rr, 1, equations)
 
             for v in eachvariable(equations)
-                cache.elements.surface_flux_values[v, 2, left_element] = f1[v]
-                cache.elements.surface_flux_values[v, 1, element] = f1[v]
+                surface_flux_values[v, 2, left_element] = f1[v]
+                surface_flux_values[v, 1, element] = f1[v]
             end
         end
     end
 
     return nothing
 end
 
-function calc_boundary_flux!(cache, u, t, boundary_conditions::NamedTuple,
+function calc_boundary_flux!(cache, t, boundary_conditions::NamedTuple,
                              mesh::StructuredMesh{1}, equations, surface_integral,
                              dg::DG)
     @unpack surface_flux = surface_integral
-    @unpack surface_flux_values, node_coordinates = cache.elements
+    @unpack surface_flux_values, node_coordinates, interfaces_u = cache.elements
+    # Boundary values are for `StructuredMesh` stored in the interface datastructure
+    boundaries_u = interfaces_u
 
     orientation = 1
 
     # Negative x-direction
     direction = 1
 
-    u_rr = get_node_vars(u, equations, dg, 1, 1)
+    u_rr = get_node_vars(boundaries_u, equations, dg, direction, 1)
     x = get_node_coords(node_coordinates, equations, dg, 1, 1)
 
     flux = boundary_conditions[direction](u_rr, orientation, direction, x, t,
                                           surface_flux, equations)
 
+    # Only flux contribution for left element, left face is the boundary flux
     for v in eachvariable(equations)
         surface_flux_values[v, direction, 1] = flux[v]
     end
 
     # Positive x-direction
     direction = 2
 
-    u_rr = get_node_vars(u, equations, dg, nnodes(dg), nelements(dg, cache))
+    u_rr = get_node_vars(boundaries_u, equations, dg, direction, nelements(dg, cache))
     x = get_node_coords(node_coordinates, equations, dg, nnodes(dg),
                         nelements(dg, cache))
 
     flux = boundary_conditions[direction](u_rr, orientation, direction, x, t,
                                           surface_flux, equations)
 
-    # Copy flux to left and right element storage
+    # Only flux contribution for right element, right face is the boundary flux
     for v in eachvariable(equations)
         surface_flux_values[v, direction, nelements(dg, cache)] = flux[v]
     end