diff --git a/src/auxiliary/precompile.jl b/src/auxiliary/precompile.jl index 5fec5875b97..335e1fa158a 100644 --- a/src/auxiliary/precompile.jl +++ b/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}}) @@ -438,7 +440,8 @@ function _precompile_manual_() # 3D, serial @assert Base.precompile(Tuple{typeof(Trixi.init_boundaries), Array{Int, 1}, TreeMesh{3, Trixi.SerialTree{3}, RealT}, - Trixi.TreeElementContainer3D{RealT, uEltype}}) + Trixi.TreeElementContainer3D{RealT, uEltype}, + basis_type_dgsem(RealT, nnodes_)}) @assert Base.precompile(Tuple{typeof(Trixi.init_interfaces), Array{Int, 1}, TreeMesh{3, Trixi.SerialTree{3}, RealT}, Trixi.TreeElementContainer3D{RealT, uEltype}}) diff --git a/src/solvers/dgsem_tree/containers_1d.jl b/src/solvers/dgsem_tree/containers_1d.jl index 0385a10ce84..b699f3ba13c 100644 --- a/src/solvers/dgsem_tree/containers_1d.jl +++ b/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 diff --git a/src/solvers/dgsem_tree/containers_2d.jl b/src/solvers/dgsem_tree/containers_2d.jl index 4e215cfd333..c7d7b71b432 100644 --- a/src/solvers/dgsem_tree/containers_2d.jl +++ b/src/solvers/dgsem_tree/containers_2d.jl @@ -353,7 +353,7 @@ 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, @@ -361,7 +361,7 @@ function init_boundaries(cell_ids, mesh::TreeMesh2D, 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 diff --git a/src/solvers/dgsem_tree/containers_3d.jl b/src/solvers/dgsem_tree/containers_3d.jl index 1b6f5de4118..da398ef1938 100644 --- a/src/solvers/dgsem_tree/containers_3d.jl +++ b/src/solvers/dgsem_tree/containers_3d.jl @@ -350,7 +350,7 @@ end # Create boundaries container and initialize boundary data in `elements`. function init_boundaries(cell_ids, mesh::TreeMesh3D, - elements::TreeElementContainer3D) + elements::TreeElementContainer3D, basis) # Initialize container n_boundaries = count_required_boundaries(mesh, cell_ids) boundaries = TreeBoundaryContainer3D{real(elements), eltype(elements)}(n_boundaries, @@ -358,7 +358,7 @@ function init_boundaries(cell_ids, mesh::TreeMesh3D, nnodes(elements)) # Connect elements with boundaries - init_boundaries!(boundaries, elements, mesh) + init_boundaries!(boundaries, elements, mesh, basis) return boundaries end @@ -388,7 +388,7 @@ function count_required_boundaries(mesh::TreeMesh3D, cell_ids) end # Initialize connectivity between elements and boundaries -function init_boundaries!(boundaries, elements, mesh::TreeMesh3D) +function init_boundaries!(boundaries, elements, mesh::TreeMesh3D, basis) # Reset boundaries count count = 0 diff --git a/src/solvers/dgsem_tree/dg_1d.jl b/src/solvers/dgsem_tree/dg_1d.jl index b4c231295b2..750b9b444ea 100644 --- a/src/solvers/dgsem_tree/dg_1d.jl +++ b/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]) diff --git a/src/solvers/dgsem_tree/dg_2d.jl b/src/solvers/dgsem_tree/dg_2d.jl index 7d2efbb7238..9e0e2ad0c36 100644 --- a/src/solvers/dgsem_tree/dg_2d.jl +++ b/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 diff --git a/src/solvers/dgsem_tree/dg_2d_parallel.jl b/src/solvers/dgsem_tree/dg_2d_parallel.jl index 37ba51e15bc..d8c23924791 100644 --- a/src/solvers/dgsem_tree/dg_2d_parallel.jl +++ b/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) diff --git a/test/test_tree_2d_acoustics.jl b/test/test_tree_2d_acoustics.jl index e38ed62cc07..49e183fdc7a 100644 --- a/test/test_tree_2d_acoustics.jl +++ b/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, diff --git a/test/test_tree_2d_euler.jl b/test/test_tree_2d_euler.jl index 5fe145d58c0..048314abff5 100644 --- a/test/test_tree_2d_euler.jl +++ b/test/test_tree_2d_euler.jl @@ -242,6 +242,29 @@ end @test_allocations(Trixi.rhs!, semi, sol, 1000) end +@trixi_testset "elixir_euler_source_terms_nonperiodic.jl" 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), + cfl=0.8, + l2=[ + 8.565448573947783e-7, + 9.279921990156959e-7, + 9.279921990210634e-7, + 2.6853435359565158e-6 + ], + linf=[ + 3.699190303185773e-6, + 4.467127135754367e-6, + 4.4671271295371184e-6, + 1.5194716922017903e-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_euler_ec.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_ec.jl"), l2=[