diff --git a/examples/tree_2d_dgsem/elixir_advection_implicit_sparse_jacobian.jl b/examples/tree_2d_dgsem/elixir_advection_implicit_sparse_jacobian.jl index abf68d33908..710bcfc66df 100644 --- a/examples/tree_2d_dgsem/elixir_advection_implicit_sparse_jacobian.jl +++ b/examples/tree_2d_dgsem/elixir_advection_implicit_sparse_jacobian.jl @@ -1,7 +1,8 @@ using Trixi using SparseConnectivityTracer # For obtaining the Jacobian sparsity pattern using SparseMatrixColorings # For obtaining the coloring vector -using OrdinaryDiffEqSDIRK, ADTypes +using OrdinaryDiffEqSDIRK, OrdinaryDiffEqDifferentiation +using ADTypes ############################################################################### ### equation, solver, mesh ### diff --git a/src/solvers/dgsem_structured/containers.jl b/src/solvers/dgsem_structured/containers.jl index fc0ad1e4776..dd6adfbee5b 100644 --- a/src/solvers/dgsem_structured/containers.jl +++ b/src/solvers/dgsem_structured/containers.jl @@ -9,6 +9,10 @@ struct StructuredElementContainer{NDIMS, RealT <: Real, uEltype <: Real, NDIMSP1, NDIMSP2, NDIMSP3} <: AbstractElementContainer # Physical coordinates at each node node_coordinates::Array{RealT, NDIMSP2} # [orientation, node_i, node_j, node_k, element] + + # Physical coordinates at boundary nodes + boundary_node_coordinates::Array{RealT, NDIMSP1} # [orientation, node_i, node_j, direction/face] + # ID of neighbor element in negative direction in orientation left_neighbors::Array{Int, 2} # [orientation, elements] @@ -39,6 +43,10 @@ function init_elements(mesh::Union{StructuredMesh{NDIMS, RealT}, node_coordinates = Array{RealT, NDIMS + 2}(undef, NDIMS, ntuple(_ -> nnodes(basis), NDIMS)..., nelements) + boundary_node_coordinates = Array{RealT, NDIMS + 1}(undef, NDIMS, + ntuple(_ -> nnodes(basis), + NDIMS - 1)..., + NDIMS * 2) left_neighbors = Array{Int, 2}(undef, NDIMS, nelements) jacobian_matrix = Array{RealT, NDIMS + 3}(undef, NDIMS, NDIMS, ntuple(_ -> nnodes(basis), NDIMS)..., @@ -58,6 +66,7 @@ function init_elements(mesh::Union{StructuredMesh{NDIMS, RealT}, elements = StructuredElementContainer{NDIMS, RealT, uEltype, NDIMS + 1, NDIMS + 2, NDIMS + 3}(node_coordinates, + boundary_node_coordinates, left_neighbors, jacobian_matrix, contravariant_vectors, diff --git a/src/solvers/dgsem_structured/containers_1d.jl b/src/solvers/dgsem_structured/containers_1d.jl index 284fdef593e..42922059c97 100644 --- a/src/solvers/dgsem_structured/containers_1d.jl +++ b/src/solvers/dgsem_structured/containers_1d.jl @@ -7,7 +7,7 @@ # Initialize data structures in element container function init_elements!(elements, mesh::StructuredMesh{1}, basis::AbstractBasisSBP) - @unpack node_coordinates, left_neighbors, + @unpack node_coordinates, boundary_node_coordinates, left_neighbors, jacobian_matrix, contravariant_vectors, inverse_jacobian = elements # Calculate node coordinates, Jacobian matrix, and inverse Jacobian determinant @@ -23,6 +23,39 @@ function init_elements!(elements, mesh::StructuredMesh{1}, basis::AbstractBasisS fill!(contravariant_vectors, NaN) initialize_left_neighbor_connectivity!(left_neighbors, mesh) + calc_boundary_node_coordinates!(boundary_node_coordinates, node_coordinates, + mesh, basis) + + return nothing +end + +function calc_boundary_node_coordinates!(boundary_node_coordinates, + node_coordinates, + mesh::StructuredMesh{1}, + basis::LobattoLegendreBasis) + nelements = size(mesh, 1) + + dim = 1 # spatial dimension + boundary_node_coordinates[dim, 1] = node_coordinates[dim, 1, 1] + boundary_node_coordinates[dim, 2] = node_coordinates[dim, nnodes(basis), nelements] + + return nothing +end + +function calc_boundary_node_coordinates!(boundary_node_coordinates, + node_coordinates, + mesh::StructuredMesh{1}, + basis::GaussLegendreBasis) + nelements = size(mesh, 1) + boundary_matrix = basis.boundary_interpolation + + dim = 1 # spatial dimension + # For structured mesh: + # Left/right boundaries are really left(-1)/right(+1) [first/second column of boundary matrix] + @views boundary_node_coordinates[dim, 1] = dot(boundary_matrix[:, 1], + node_coordinates[dim, :, 1]) + @views boundary_node_coordinates[dim, 2] = dot(boundary_matrix[:, 2], + node_coordinates[dim, :, nelements]) return nothing end diff --git a/src/solvers/dgsem_structured/dg_1d.jl b/src/solvers/dgsem_structured/dg_1d.jl index 781c39dfcbf..6c3266137da 100644 --- a/src/solvers/dgsem_structured/dg_1d.jl +++ b/src/solvers/dgsem_structured/dg_1d.jl @@ -78,7 +78,7 @@ 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, interfaces_u = cache.elements + @unpack surface_flux_values, boundary_node_coordinates, interfaces_u = cache.elements # Boundary values are for `StructuredMesh` stored in the interface datastructure boundaries_u = interfaces_u @@ -88,7 +88,7 @@ function calc_boundary_flux!(cache, t, boundary_conditions::NamedTuple, direction = 1 u_rr = get_node_vars(boundaries_u, equations, dg, direction, 1) - x = get_node_coords(node_coordinates, equations, dg, 1, 1) + x = get_node_coords(boundary_node_coordinates, equations, dg, direction) flux = boundary_conditions[direction](u_rr, orientation, direction, x, t, surface_flux, equations) @@ -102,8 +102,7 @@ function calc_boundary_flux!(cache, t, boundary_conditions::NamedTuple, direction = 2 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)) + x = get_node_coords(boundary_node_coordinates, equations, dg, direction) flux = boundary_conditions[direction](u_rr, orientation, direction, x, t, surface_flux, equations) diff --git a/src/solvers/dgsem_tree/containers_1d.jl b/src/solvers/dgsem_tree/containers_1d.jl index b699f3ba13c..15ddff369be 100644 --- a/src/solvers/dgsem_tree/containers_1d.jl +++ b/src/solvers/dgsem_tree/containers_1d.jl @@ -374,9 +374,9 @@ function count_required_boundaries(mesh::TreeMesh1D, cell_ids) 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) +function calc_boundary_node_coordinates!(boundaries, element, count, direction, + elements, mesh::TreeMesh1D, + basis::LobattoLegendreBasis) el_node_coords = elements.node_coordinates bnd_node_coords = boundaries.node_coordinates @@ -395,9 +395,9 @@ function set_boundary_node_coordinates!(boundaries, element, count, direction, 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) +function calc_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 @@ -465,9 +465,9 @@ function init_boundaries!(boundaries, elements, mesh::TreeMesh1D, basis) # Set orientation (x -> 1) boundaries.orientations[count] = 1 - # Store node coordinates - set_boundary_node_coordinates!(boundaries, element, count, direction, - elements, mesh, basis) + # Calculate node coordinates + calc_boundary_node_coordinates!(boundaries, element, count, direction, + elements, mesh, basis) end end diff --git a/src/solvers/dgsem_tree/containers_2d.jl b/src/solvers/dgsem_tree/containers_2d.jl index 6c5d970f70c..be76b91c102 100644 --- a/src/solvers/dgsem_tree/containers_2d.jl +++ b/src/solvers/dgsem_tree/containers_2d.jl @@ -391,9 +391,9 @@ function count_required_boundaries(mesh::TreeMesh2D, cell_ids) 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) +function calc_boundary_node_coordinates!(boundaries, element, count, direction, + elements, mesh::TreeMesh2D, + basis::LobattoLegendreBasis) el_node_coords = elements.node_coordinates bnd_node_coords = boundaries.node_coordinates @@ -413,9 +413,9 @@ function set_boundary_node_coordinates!(boundaries, element, count, direction, 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) +function calc_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 @@ -527,9 +527,9 @@ function init_boundaries!(boundaries, elements, mesh::TreeMesh2D, basis) boundaries.orientations[count] = 2 # y direction end - # Store node coordinates - set_boundary_node_coordinates!(boundaries, element, count, direction, - elements, mesh, basis) + # Calculate node coordinates + calc_boundary_node_coordinates!(boundaries, element, count, direction, + elements, mesh, basis) end end diff --git a/test/test_structured_1d.jl b/test/test_structured_1d.jl index cacc799b148..74f12927117 100644 --- a/test/test_structured_1d.jl +++ b/test/test_structured_1d.jl @@ -201,6 +201,27 @@ 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), + # Identical errors as for the `TreeMesh` version of this example + 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_source_terms_nonperiodic_fvO2.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms_nonperiodic_fvO2.jl"),