Gauss-Legendre: Boundary Coordinates for 1D TreeMesh

src/auxiliary/precompile.jl

@@ -396,7 +396,8 @@ function _precompile_manual_()
         # 1D, serial
         @assert Base.precompile(Tuple{typeof(Trixi.init_boundaries), Array{Int, 1},
                                       TreeMesh{1, Trixi.SerialTree{1}, RealT},
-                                      Trixi.TreeElementContainer1D{RealT, uEltype}})
+                                      Trixi.TreeElementContainer1D{RealT, uEltype},
+                                      basis_type_dgsem(RealT, nnodes_)})
         @assert Base.precompile(Tuple{typeof(Trixi.init_interfaces), Array{Int, 1},
                                       TreeMesh{1, Trixi.SerialTree{1}, RealT},
                                       Trixi.TreeElementContainer1D{RealT, uEltype}})

src/solvers/dgsem_tree/containers_1d.jl

@@ -337,14 +337,14 @@ end
 
 # Create boundaries container and initialize boundary data in `elements`.
 function init_boundaries(cell_ids, mesh::TreeMesh1D,
-                         elements::TreeElementContainer1D)
+                         elements::TreeElementContainer1D, basis)
     # Initialize container
     n_boundaries = count_required_boundaries(mesh, cell_ids)
     boundaries = TreeBoundaryContainer1D{real(elements), eltype(elements)}(n_boundaries,
                                                                            nvariables(elements))
 
     # Connect elements with boundaries
-    init_boundaries!(boundaries, elements, mesh)
+    init_boundaries!(boundaries, elements, mesh, basis)
     return boundaries
 end
 
@@ -373,8 +373,50 @@ function count_required_boundaries(mesh::TreeMesh1D, cell_ids)
     return count
 end
 
+# For Lobtto points, we can simply use the outer nodes of the elements as boundary nodes.
+function set_boundary_node_coordinates!(boundaries, element, count, direction,
+                                        elements, mesh::TreeMesh1D,
+                                        basis::LobattoLegendreBasis)
+    orientation = 1 # always 1 in 1D
+    if direction == 1
+        @views boundaries.node_coordinates[:, count] .= elements.node_coordinates[orientation,
+                                                                                  1,
+                                                                                  element]
+    elseif direction == 2
+        @views boundaries.node_coordinates[:, count] .= elements.node_coordinates[orientation,
+                                                                                  end,
+                                                                                  element]
+    else
+        error("should not happen")
+    end
+
+    return nothing
+end
+
+# For Gauss points, we need to interpolate the boundary node coordinates.
+function set_boundary_node_coordinates!(boundaries, element, count, direction,
+                                        elements, mesh::TreeMesh1D,
+                                        basis::GaussLegendreBasis)
+    orientation = 1 # always 1 in 1D
+    if direction == 1
+        @views x_interpolated_left = dot(basis.boundary_interpolation[:, 1],
+                                         elements.node_coordinates[orientation, :,
+                                                                   element])
+        boundaries.node_coordinates[:, count] .= x_interpolated_left
+    elseif direction == 2
+        @views x_interpolated_right = dot(basis.boundary_interpolation[:, 2],
+                                          elements.node_coordinates[orientation, :,
+                                                                    element])
+        boundaries.node_coordinates[:, count] .= x_interpolated_right
+    else
+        error("should not happen")
+    end
+
+    return nothing
+end
+
 # Initialize connectivity between elements and boundaries
-function init_boundaries!(boundaries, elements, mesh::TreeMesh1D)
+function init_boundaries!(boundaries, elements, mesh::TreeMesh1D, basis)
     # Reset boundaries count
     count = 0
 
@@ -419,14 +461,8 @@ function init_boundaries!(boundaries, elements, mesh::TreeMesh1D)
             boundaries.orientations[count] = 1
 
             # Store node coordinates
-            enc = elements.node_coordinates
-            if direction == 1 # -x direction
-                boundaries.node_coordinates[:, count] .= enc[:, 1, element]
-            elseif direction == 2 # +x direction
-                boundaries.node_coordinates[:, count] .= enc[:, end, element]
-            else
-                error("should not happen")
-            end
+            set_boundary_node_coordinates!(boundaries, element, count, direction,
+                                           elements, mesh, basis)
         end
     end
 
@@ -462,7 +498,7 @@ function reinitialize_containers!(mesh::TreeMesh{1}, equations, dg::DGSEM, cache
     # re-initialize boundaries container
     @unpack boundaries = cache
     resize!(boundaries, count_required_boundaries(mesh, leaf_cell_ids))
-    init_boundaries!(boundaries, elements, mesh)
+    init_boundaries!(boundaries, elements, mesh, dg.basis)
 
     return nothing
 end

src/solvers/dgsem_tree/dg_1d.jl

@@ -20,7 +20,7 @@ function create_cache(mesh::TreeMesh{1}, equations,
 
     interfaces = init_interfaces(leaf_cell_ids, mesh, elements)
 
-    boundaries = init_boundaries(leaf_cell_ids, mesh, elements)
+    boundaries = init_boundaries(leaf_cell_ids, mesh, elements, dg.basis)
 
     # Container cache
     cache = (; elements, interfaces, boundaries)
@@ -417,9 +417,11 @@ function prolong2interfaces!(cache, u_or_flux_viscous,
             for ii in eachnode(dg)
                 # Not += to allow `@muladd` to turn these into FMAs
                 # (see comment at the top of the file)
+                # Need `boundary_interpolation` at right (+1) node for left element
                 interface_u_1 = (interface_u_1 +
                                  u_or_flux_viscous[v, ii, left_element] *
                                  boundary_interpolation[ii, 2])
+                # Need `boundary_interpolation` at left (-1) node for right element
                 interface_u_2 = (interface_u_2 +
                                  u_or_flux_viscous[v, ii, right_element] *
                                  boundary_interpolation[ii, 1])
@@ -531,6 +533,52 @@ function prolong2boundaries!(cache, u_or_flux_viscous,
     return nothing
 end
 
+function prolong2boundaries!(cache, u_or_flux_viscous,
+                             mesh::TreeMesh{1}, equations,
+                             dg::DG{<:GaussLegendreBasis})
+    @unpack boundaries = cache
+    @unpack neighbor_sides = boundaries
+    @unpack boundary_interpolation = dg.basis
+
+    @threaded for boundary in eachboundary(dg, cache)
+        element = boundaries.neighbor_ids[boundary]
+
+        # boundary in x-direction
+        if neighbor_sides[boundary] == 1
+            # element in -x direction of boundary => need to evaluate at right boundary node (+1)
+            for v in eachvariable(equations)
+                # Interpolate to the boundaries using a local variable for
+                # the accumulation of values (to reduce global memory operations).
+                boundary_u_1 = zero(eltype(boundaries.u))
+                for ii in eachnode(dg)
+                    # Not += to allow `@muladd` to turn these into FMAs
+                    # (see comment at the top of the file)
+                    boundary_u_1 = (boundary_u_1 +
+                                    u_or_flux_viscous[v, ii, element] *
+                                    boundary_interpolation[ii, 2])
+                end
+                boundaries.u[1, v, boundary] = boundary_u_1
+            end
+        else # Element in +x direction of boundary => need to evaluate at left boundary node (-1)
+            for v in eachvariable(equations)
+                # Interpolate to the boundaries using a local variable for
+                # the accumulation of values (to reduce global memory operations).
+                boundary_u_2 = zero(eltype(boundaries.u))
+                for ii in eachnode(dg)
+                    # Not += to allow `@muladd` to turn these into FMAs
+                    # (see comment at the top of the file)
+                    boundary_u_2 = (boundary_u_2 +
+                                    u_or_flux_viscous[v, ii, element] *
+                                    boundary_interpolation[ii, 1])
+                end
+                boundaries.u[2, v, boundary] = boundary_u_2
+            end
+        end
+    end
+
+    return nothing
+end
+
 function calc_boundary_flux!(cache, t, boundary_conditions::NamedTuple,
                              mesh::TreeMesh{1}, equations, surface_integral, dg::DG)
     @unpack surface_flux_values = cache.elements

test/test_tree_1d_euler.jl

@@ -147,6 +147,26 @@ end
     @test_allocations(Trixi.rhs!, semi, sol, 1000)
 end
 
+@trixi_testset "elixir_euler_source_terms_nonperiodic.jl (Gauss-Legendre)" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR,
+                                 "elixir_euler_source_terms_nonperiodic.jl"),
+                        solver=DGSEM(polydeg = 3, basis_type = GaussLegendreBasis,
+                                     surface_flux = flux_lax_friedrichs),
+                        l2=[
+                            6.179119971404758e-7,
+                            6.831335637140733e-7,
+                            1.8153512648336213e-6
+                        ],
+                        linf=[
+                            2.3035825069683824e-6,
+                            2.7398314812465685e-6,
+                            7.132056524916663e-6
+                        ])
+    # 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_euler_ec.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_ec.jl"),
                         l2=[