diff --git a/src/auxiliary/precompile.jl b/src/auxiliary/precompile.jl index 7f5cec6e731..5fec5875b97 100644 --- a/src/auxiliary/precompile.jl +++ b/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}}) diff --git a/src/solvers/dgsem_tree/containers_1d.jl b/src/solvers/dgsem_tree/containers_1d.jl index 0f562670756..0385a10ce84 100644 --- a/src/solvers/dgsem_tree/containers_1d.jl +++ b/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 diff --git a/src/solvers/dgsem_tree/dg_1d.jl b/src/solvers/dgsem_tree/dg_1d.jl index bfbd3784bd0..b4c231295b2 100644 --- a/src/solvers/dgsem_tree/dg_1d.jl +++ b/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 diff --git a/test/test_tree_1d_euler.jl b/test/test_tree_1d_euler.jl index 42889430bf0..c786afb5a52 100644 --- a/test/test_tree_1d_euler.jl +++ b/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=[