Gauss-Legendre DGSEM Solver: Structured 1D: Boundaries

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 ###

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,

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

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)

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
 

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
 

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"),