Gauss-Legendre DGSEM Solver: Boundary Nodes for 2D TreeMesh

src/auxiliary/precompile.jl

@@ -407,7 +407,8 @@ function _precompile_manual_()
         # 2D, serial
         @assert Base.precompile(Tuple{typeof(Trixi.init_boundaries), Array{Int, 1},
                                       TreeMesh{2, Trixi.SerialTree{2}, RealT},
-                                      Trixi.TreeElementContainer2D{RealT, uEltype}})
+                                      Trixi.TreeElementContainer2D{RealT, uEltype},
+                                      basis_type_dgsem(RealT, nnodes_)})
         @assert Base.precompile(Tuple{typeof(Trixi.init_interfaces), Array{Int, 1},
                                       TreeMesh{2, Trixi.SerialTree{2}, RealT},
                                       Trixi.TreeElementContainer2D{RealT, uEltype}})
@@ -421,7 +422,8 @@ function _precompile_manual_()
         # 2D, parallel
         @assert Base.precompile(Tuple{typeof(Trixi.init_boundaries), Array{Int, 1},
                                       TreeMesh{2, Trixi.ParallelTree{2}, RealT},
-                                      Trixi.TreeElementContainer2D{RealT, uEltype}})
+                                      Trixi.TreeElementContainer2D{RealT, uEltype},
+                                      basis_type_dgsem(RealT, nnodes_)})
         @assert Base.precompile(Tuple{typeof(Trixi.init_interfaces), Array{Int, 1},
                                       TreeMesh{2, Trixi.ParallelTree{2}, RealT},
                                       Trixi.TreeElementContainer2D{RealT, uEltype}})

src/solvers/dgsem_tree/containers_1d.jl

@@ -377,15 +377,16 @@ end
 function set_boundary_node_coordinates!(boundaries, element, count, direction,
                                         elements, mesh::TreeMesh1D,
                                         basis::LobattoLegendreBasis)
+    el_node_coords = elements.node_coordinates
+    bnd_node_coords = boundaries.node_coordinates
+
     orientation = 1 # always 1 in 1D
     if direction == 1
-        @views boundaries.node_coordinates[:, count] .= elements.node_coordinates[orientation,
-                                                                                  1,
-                                                                                  element]
+        bnd_node_coords[orientation, count] = el_node_coords[orientation, 1,
+                                                             element]
     elseif direction == 2
-        @views boundaries.node_coordinates[:, count] .= elements.node_coordinates[orientation,
-                                                                                  end,
-                                                                                  element]
+        bnd_node_coords[orientation, count] = el_node_coords[orientation, end,
+                                                             element]
     else
         error("should not happen")
     end
@@ -397,17 +398,21 @@ end
 function set_boundary_node_coordinates!(boundaries, element, count, direction,
                                         elements, mesh::TreeMesh1D,
                                         basis::GaussLegendreBasis)
+    boundary_matrix = basis.boundary_interpolation
+    el_node_coords = elements.node_coordinates
+    bnd_node_coords = boundaries.node_coordinates
+
     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
+        @views x_interpolated_left = dot(boundary_matrix[:, 1],
+                                         el_node_coords[orientation, :,
+                                                        element])
+        bnd_node_coords[orientation, 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
+        @views x_interpolated_right = dot(boundary_matrix[:, 2],
+                                          el_node_coords[orientation, :,
+                                                         element])
+        bnd_node_coords[orientation, count] = x_interpolated_right
     else
         error("should not happen")
     end

src/solvers/dgsem_tree/containers_2d.jl

@@ -353,15 +353,15 @@ end
 
 # Create boundaries container and initialize boundary data in `elements`.
 function init_boundaries(cell_ids, mesh::TreeMesh2D,
-                         elements::TreeElementContainer2D)
+                         elements::TreeElementContainer2D, basis)
     # Initialize container
     n_boundaries = count_required_boundaries(mesh, cell_ids)
     boundaries = TreeBoundaryContainer2D{real(elements), eltype(elements)}(n_boundaries,
                                                                            nvariables(elements),
                                                                            nnodes(elements))
 
     # Connect elements with boundaries
-    init_boundaries!(boundaries, elements, mesh)
+    init_boundaries!(boundaries, elements, mesh, basis)
     return boundaries
 end
 
@@ -390,8 +390,89 @@ function count_required_boundaries(mesh::TreeMesh2D, cell_ids)
     return count
 end
 
+# For Lobatto 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::TreeMesh2D,
+                                        basis::LobattoLegendreBasis)
+    el_node_coords = elements.node_coordinates
+    bnd_node_coords = boundaries.node_coordinates
+
+    if direction == 1 # -x direction
+        @views bnd_node_coords[:, :, count] .= el_node_coords[:, 1, :, element]
+    elseif direction == 2 # +x direction
+        @views bnd_node_coords[:, :, count] .= el_node_coords[:, end, :, element]
+    elseif direction == 3 # -y direction
+        @views bnd_node_coords[:, :, count] .= el_node_coords[:, :, 1, element]
+    elseif direction == 4 # +y direction
+        @views bnd_node_coords[:, :, count] .= el_node_coords[:, :, 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::TreeMesh2D,
+                                        basis::GaussLegendreBasis)
+    boundary_matrix = basis.boundary_interpolation
+    el_node_coords = elements.node_coordinates
+    bnd_node_coords = boundaries.node_coordinates
+
+    if direction == 1 # -x direction: interpolate in x for each y node j
+        for j in eachnode(basis)
+            for orientation in 1:2 # Need to set both x and y coordinate of boundary node
+                @views bnd_node_coords[orientation, j, count] = dot(boundary_matrix[:,
+                                                                                    1],
+                                                                    el_node_coords[orientation,
+                                                                                   :,
+                                                                                   j,
+                                                                                   element])
+            end
+        end
+    elseif direction == 2 # +x direction: interpolate in x for each y node j
+        for j in eachnode(basis)
+            for orientation in 1:2 # Need to set both x and y coordinate of boundary node
+                @views bnd_node_coords[orientation, j, count] = dot(boundary_matrix[:,
+                                                                                    2],
+                                                                    el_node_coords[orientation,
+                                                                                   :,
+                                                                                   j,
+                                                                                   element])
+            end
+        end
+    elseif direction == 3 # -y direction: interpolate in y for each x node i
+        for i in eachnode(basis)
+            for orientation in 1:2 # Need to set both x and y coordinate of boundary node
+                @views bnd_node_coords[orientation, i, count] = dot(boundary_matrix[:,
+                                                                                    1],
+                                                                    el_node_coords[orientation,
+                                                                                   i,
+                                                                                   :,
+                                                                                   element])
+            end
+        end
+    elseif direction == 4 # +y direction: interpolate in y for each x node i
+        for i in eachnode(basis)
+            for orientation in 1:2 # Need to set both x and y coordinate of boundary node
+                @views bnd_node_coords[orientation, i, count] = dot(boundary_matrix[:,
+                                                                                    2],
+                                                                    el_node_coords[orientation,
+                                                                                   i,
+                                                                                   :,
+                                                                                   element])
+            end
+        end
+    else
+        error("should not happen")
+    end
+
+    return nothing
+end
+
 # Initialize connectivity between elements and boundaries
-function init_boundaries!(boundaries, elements, mesh::TreeMesh2D)
+function init_boundaries!(boundaries, elements, mesh::TreeMesh2D, basis)
     # Exit early if there are no boundaries to initialize
     if nboundaries(boundaries) == 0
         # In this case n_boundaries_per_direction still needs to be reset!
@@ -441,24 +522,14 @@ function init_boundaries!(boundaries, elements, mesh::TreeMesh2D)
 
             # Set orientation (x -> 1, y -> 2)
             if direction in (1, 2)
-                boundaries.orientations[count] = 1
+                boundaries.orientations[count] = 1 # x direction
             else
-                boundaries.orientations[count] = 2
+                boundaries.orientations[count] = 2 # y direction
             end
 
             # 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]
-            elseif direction == 3 # -y direction
-                boundaries.node_coordinates[:, :, count] .= enc[:, :, 1, element]
-            elseif direction == 4 # +y 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
 
@@ -1349,7 +1420,7 @@ function reinitialize_containers!(mesh::Union{TreeMesh{2}, TreeMesh{3}}, equatio
     # 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)
 
     # re-initialize mortars container
     @unpack mortars = cache

src/solvers/dgsem_tree/dg_1d.jl

@@ -561,12 +561,8 @@ function prolong2boundaries!(cache, u_or_flux_viscous,
             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])

src/solvers/dgsem_tree/dg_2d.jl

@@ -20,7 +20,7 @@ function create_cache(mesh::Union{TreeMesh{2}, TreeMesh{3}}, 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)
 
     mortars = init_mortars(leaf_cell_ids, mesh, elements, dg.mortar)
 
@@ -699,6 +699,74 @@ function prolong2boundaries!(cache, u,
     return nothing
 end
 
+function prolong2boundaries!(cache, u,
+                             mesh::TreeMesh{2}, equations,
+                             dg::DG{<:GaussLegendreBasis})
+    @unpack boundaries = cache
+    @unpack orientations, neighbor_sides = boundaries
+    @unpack boundary_interpolation = dg.basis
+
+    @threaded for boundary in eachboundary(dg, cache)
+        element = boundaries.neighbor_ids[boundary]
+
+        if orientations[boundary] == 1
+            # boundary in x-direction
+            if neighbor_sides[boundary] == 1
+                # element in -x direction of boundary => interpolate to right boundary node (+1)
+                for l in eachnode(dg), v in eachvariable(equations)
+                    # Interpolate to the boundaries using a local variable for
+                    # the accumulation of values (to reduce global memory operations).
+                    boundary_u = 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 = (boundary_u +
+                                      u[v, ii, l, element] *
+                                      boundary_interpolation[ii, 2])
+                    end
+                    boundaries.u[1, v, l, boundary] = boundary_u
+                end
+            else # element in +x direction of boundary => interpolate to left boundary node (-1)
+                for l in eachnode(dg), v in eachvariable(equations)
+                    boundary_u = zero(eltype(boundaries.u))
+                    for ii in eachnode(dg)
+                        boundary_u = (boundary_u +
+                                      u[v, ii, l, element] *
+                                      boundary_interpolation[ii, 1])
+                    end
+                    boundaries.u[2, v, l, boundary] = boundary_u
+                end
+            end
+        else # if orientations[boundary] == 2
+            # boundary in y-direction
+            if neighbor_sides[boundary] == 1
+                # element in -y direction of boundary => interpolate to right boundary node (+1)
+                for l in eachnode(dg), v in eachvariable(equations)
+                    boundary_u = zero(eltype(boundaries.u))
+                    for ii in eachnode(dg)
+                        boundary_u = (boundary_u +
+                                      u[v, l, ii, element] *
+                                      boundary_interpolation[ii, 2])
+                    end
+                    boundaries.u[1, v, l, boundary] = boundary_u
+                end
+            else # element in +y direction of boundary => interpolate to left boundary node (-1)
+                for l in eachnode(dg), v in eachvariable(equations)
+                    boundary_u = zero(eltype(boundaries.u))
+                    for ii in eachnode(dg)
+                        boundary_u = (boundary_u +
+                                      u[v, l, ii, element] *
+                                      boundary_interpolation[ii, 1])
+                    end
+                    boundaries.u[2, v, l, boundary] = boundary_u
+                end
+            end
+        end
+    end
+
+    return nothing
+end
+
 function calc_boundary_flux!(cache, t, boundary_conditions::NamedTuple,
                              mesh::TreeMesh{2}, equations, surface_integral, dg::DG)
     @unpack surface_flux_values = cache.elements

src/solvers/dgsem_tree/dg_2d_parallel.jl

@@ -253,7 +253,7 @@ function create_cache(mesh::TreeMeshParallel{2}, equations,
 
     mpi_interfaces = init_mpi_interfaces(leaf_cell_ids, mesh, elements)
 
-    boundaries = init_boundaries(leaf_cell_ids, mesh, elements)
+    boundaries = init_boundaries(leaf_cell_ids, mesh, elements, dg.basis)
 
     mortars = init_mortars(leaf_cell_ids, mesh, elements, dg.mortar)
 

test/test_tree_2d_acoustics.jl

@@ -87,19 +87,55 @@ end
 
 @trixi_testset "elixir_acoustics_gauss_wall.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_acoustics_gauss_wall.jl"),
-                        l2=[0.019419398248465843, 0.019510701017551826,
+                        l2=[
+                            0.019419398248465843,
+                            0.019510701017551826,
                             0.04818246051887614,
-                            7.382060834820337e-17, 0.0, 1.4764121669640674e-16,
+                            7.382060834820337e-17,
+                            0.0,
+                            1.4764121669640674e-16,
                             1.4764121669640674e-16],
-                        linf=[0.18193631937316496, 0.1877464607867628,
+                        linf=[
+                            0.18193631937316496,
+                            0.1877464607867628,
                             1.0355388011792845,
-                            2.220446049250313e-16, 0.0, 4.440892098500626e-16,
+                            2.220446049250313e-16,
+                            0.0,
+                            4.440892098500626e-16,
                             4.440892098500626e-16])
     # 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_acoustics_gauss_wall.jl (Gauss Legendre)" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_acoustics_gauss_wall.jl"),
+                        solver=DGSEM(polydeg = 5, surface_flux = flux_lax_friedrichs,
+                                     basis_type = GaussLegendreBasis),
+                        cfl=0.6,
+                        l2=[
+                            0.01944153623864891,
+                            0.01952877141847981,
+                            0.04820571764883919,
+                            1.1071998298551595e-16,
+                            0.0,
+                            2.214399659710319e-16,
+                            2.214399659710319e-16
+                        ],
+                        linf=[
+                            0.1828989576562236,
+                            0.18857385551148917,
+                            1.036543390095062,
+                            3.3306690738754696e-16,
+                            0.0,
+                            6.661338147750939e-16,
+                            6.661338147750939e-16
+                        ])
+    # 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_acoustics_monopole.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_acoustics_monopole.jl"),
                         l2=[0.006816790293009947, 0.0065068948357351625,