Sort unstructured boundary conditions

src/equations/compressible_euler_2d.jl

@@ -592,7 +592,6 @@ velocity of the the exterior fictitious element to the negative of the internal
   normal = normal_direction / norm(normal_direction)
 
   # compute the normal and tangential components of the velocity
-  u_tangent = zeros(eltype(normal), 2)
   u_normal  = normal[1] * u_internal[2] + normal[2] * u_internal[3]
   u_tangent = (u_internal[2] - u_normal * normal[1], u_internal[3] - u_normal * normal[2])
 

src/mesh/mesh_io.jl

@@ -125,6 +125,7 @@ function save_mesh_file(mesh::UnstructuredQuadMesh, output_directory)
     attributes(file)["ndims"] = ndims(mesh)
     attributes(file)["size"] = length(mesh)
     attributes(file)["mesh_filename"] = mesh.filename
+    attributes(file)["periodicity"] = collect(mesh.periodicity)
   end
 
   return filename
@@ -197,12 +198,11 @@ function load_mesh_serial(mesh_file::AbstractString; n_cells_max, RealT)
 
     mesh = CurvedMesh(size, mapping; RealT=RealT, unsaved_changes=false, mapping_as_string=mapping_as_string)
   elseif mesh_type == "UnstructuredQuadMesh"
-    # FIXME: This works for visualization purposes via Trixi2Vtk however one currently cannot
-    #        restrat an unstructured and periodic simulation
-    mesh_filename = h5open(mesh_file, "r") do file
-      return read(attributes(file)["mesh_filename"])
+    mesh_filename, periodicity_ = h5open(mesh_file, "r") do file
+      return read(attributes(file)["mesh_filename"]),
+             read(attributes(file)["periodicity"])
     end
-    mesh = UnstructuredQuadMesh(mesh_filename; RealT=RealT, periodicity=false, unsaved_changes=false)
+    mesh = UnstructuredQuadMesh(mesh_filename; RealT=RealT, periodicity=periodicity_, unsaved_changes=false)
   else
     error("Unknown mesh type!")
   end

src/semidiscretization/semidiscretization_hyperbolic.jl

@@ -52,7 +52,7 @@ function SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver
                                       initial_cache=NamedTuple())
 
   cache = (; create_cache(mesh, equations, solver, RealT, uEltype)..., initial_cache...)
-  _boundary_conditions = digest_boundary_conditions(boundary_conditions)
+  _boundary_conditions = digest_boundary_conditions(boundary_conditions, cache)
 
   SemidiscretizationHyperbolic{typeof(mesh), typeof(equations), typeof(initial_condition), typeof(_boundary_conditions), typeof(source_terms), typeof(solver), typeof(cache)}(
     mesh, equations, initial_condition, _boundary_conditions, source_terms, solver, cache)
@@ -81,22 +81,27 @@ end
 
 
 # allow passing named tuples of BCs constructed in an arbitrary order
-digest_boundary_conditions(boundary_conditions) = boundary_conditions
+digest_boundary_conditions(boundary_conditions, cache) = boundary_conditions
 
-function digest_boundary_conditions(boundary_conditions::NamedTuple{Keys,ValueTypes}) where {Keys, ValueTypes<:NTuple{2,Any}} # 1D
+function digest_boundary_conditions(boundary_conditions::NamedTuple{Keys,ValueTypes}, cache) where {Keys, ValueTypes<:NTuple{2,Any}} # 1D
   @unpack x_neg, x_pos = boundary_conditions
   (; x_neg, x_pos)
 end
-function digest_boundary_conditions(boundary_conditions::NamedTuple{Keys,ValueTypes}) where {Keys, ValueTypes<:NTuple{4,Any}} # 2D
+function digest_boundary_conditions(boundary_conditions::NamedTuple{Keys,ValueTypes}, cache) where {Keys, ValueTypes<:NTuple{4,Any}} # 2D
   @unpack x_neg, x_pos, y_neg, y_pos = boundary_conditions
   (; x_neg, x_pos, y_neg, y_pos)
 end
-function digest_boundary_conditions(boundary_conditions::NamedTuple{Keys,ValueTypes}) where {Keys, ValueTypes<:NTuple{6,Any}} # 3D
+function digest_boundary_conditions(boundary_conditions::NamedTuple{Keys,ValueTypes}, cache) where {Keys, ValueTypes<:NTuple{6,Any}} # 3D
   @unpack x_neg, x_pos, y_neg, y_pos, z_neg, z_pos = boundary_conditions
   (; x_neg, x_pos, y_neg, y_pos, z_neg, z_pos)
 end
 
-function digest_boundary_conditions(boundary_conditions::AbstractArray)
+# sort the boundary conditions from a dictionary and into tuples
+function digest_boundary_conditions(boundary_conditions::Dict, cache)
+  UnstructuredQuadSortedBoundaryTypes(boundary_conditions, cache)
+end
+
+function digest_boundary_conditions(boundary_conditions::AbstractArray, cache)
   throw(ArgumentError("Please use a (named) tuple instead of an (abstract) array to supply multiple boundary conditions (to improve performance)."))
 end
 
@@ -130,9 +135,10 @@ function Base.show(io::IO, ::MIME"text/plain", semi::SemidiscretizationHyperboli
     summary_line(io, "mesh", semi.mesh)
     summary_line(io, "equations", semi.equations |> typeof |> nameof)
     summary_line(io, "initial condition", semi.initial_condition)
-    if semi.boundary_conditions isa Dict
-      summary_line(io, "boundary conditions", length(semi.boundary_conditions))
-      for (boundary_name, boundary_condition) in semi.boundary_conditions
+    if semi.boundary_conditions isa UnstructuredQuadSortedBoundaryTypes
+      @unpack boundary_dictionary = semi.boundary_conditions
+      summary_line(io, "boundary conditions", length(boundary_dictionary))
+      for (boundary_name, boundary_condition) in boundary_dictionary
         summary_line(increment_indent(io), boundary_name, typeof(boundary_condition))
       end
     else # non dictionary boundary conditions container

src/solvers/dg_unstructured_quad/dg.jl

@@ -18,4 +18,5 @@ end
 include("mappings_geometry_curved_2d.jl")
 include("mappings_geometry_straight_2d.jl")
 include("containers_2d.jl")
+include("sort_boundary_conditions.jl")
 include("dg_2d.jl")

src/solvers/dg_unstructured_quad/dg_2d.jl

@@ -223,35 +223,65 @@ end
 
 function calc_boundary_flux!(cache, t, boundary_conditions,
                              mesh::UnstructuredQuadMesh, equations, dg::DG)
-  @unpack surface_flux_values = cache.elements
-  @unpack element_id, element_side_id, name  = cache.boundaries
+  @unpack boundary_condition_types, boundary_indices = boundary_conditions
 
-  @threaded for boundary in eachboundary(dg, cache)
-    # Get the element and side IDs on the boundary element
-    element = element_id[boundary]
-    side    = element_side_id[boundary]
+  calc_boundary_flux_by_type!(cache, t, boundary_condition_types, boundary_indices,
+                              mesh, equations, dg)
+  return nothing
+end
 
-    # pull the external state function from the bounary condition dictionary
-    boundary_condition = boundary_conditions[ name[boundary] ]
 
-    # calc boundary flux on the current boundary interface
-    calc_boundary_flux!(surface_flux_values, t, boundary_condition, mesh, equations, dg, cache,
-                        side, element, boundary)
-  end
+# Iterate over tuples of boundary condition types and associated indices
+# in a type-stable way using "lispy tuple programming".
+function calc_boundary_flux_by_type!(cache, t, BCs::NTuple{N,Any},
+                                     BC_indices::NTuple{N,Vector{Int}},
+                                     mesh::UnstructuredQuadMesh, equations, dg::DG) where {N}
+  # Extract the boundary condition type and index vector
+  boundary_condition = first(BCs)
+  boundary_condition_indices = first(BC_indices)
+  # Extract the remaining types and indices to be processed later
+  remaining_boundary_conditions = Base.tail(BCs)
+  remaining_boundary_condition_indices = Base.tail(BC_indices)
+
+  # process the first boundary condition type
+  calc_boundary_flux!(cache, t, boundary_condition, boundary_condition_indices, mesh, equations, dg)
+
+  # recursively call this method with the unprocessed boundary types
+  calc_boundary_flux_by_type!(cache, t, remaining_boundary_conditions,
+                              remaining_boundary_condition_indices, mesh, equations, dg)
 
   return nothing
 end
 
-# use a function barrier for now to improve type stability
-@noinline function calc_boundary_flux!(surface_flux_values, t, boundary_condition::BC,
-                                       mesh::UnstructuredQuadMesh, equations, dg::DG, cache,
-                                       side, element, boundary) where {BC}
-  for node in eachnode(dg)
-    calc_boundary_flux!(surface_flux_values, t, boundary_condition, mesh, equations, dg, cache,
-                        node, side, element, boundary)
+# terminate the type-stable iteration over tuples
+function calc_boundary_flux_by_type!(cache, t, BCs::Tuple{}, BC_indices::Tuple{},
+                                     mesh::UnstructuredQuadMesh, equations, dg::DG)
+  nothing
+end
+
+
+function calc_boundary_flux!(cache, t, boundary_condition, boundary_indexing,
+                             mesh::UnstructuredQuadMesh, equations, dg::DG)
+  @unpack surface_flux_values = cache.elements
+  @unpack element_id, element_side_id = cache.boundaries
+
+  @threaded for local_index in eachindex(boundary_indexing)
+    # use the local index to get the global boundary index from the pre-sorted list
+    boundary = boundary_indexing[local_index]
+
+    # get the element and side IDs on the boundary element
+    element = element_id[boundary]
+    side    = element_side_id[boundary]
+
+    # calc boundary flux on the current boundary interface
+    for node in eachnode(dg)
+      calc_boundary_flux!(surface_flux_values, t, boundary_condition, mesh, equations, dg, cache,
+                          node, side, element, boundary)
+    end
   end
 end
 
+
 # inlined version of the boundary flux calculation along a physical interface where the
 # boundary flux values are set according to a particular `boundary_condition` function
 @inline function calc_boundary_flux!(surface_flux_values, t, boundary_condition,

src/solvers/dg_unstructured_quad/sort_boundary_conditions.jl

@@ -0,0 +1,50 @@
+"""
+    UnstructuredQuadSortedBoundaryTypes{N, BCs<:NTuple{N, Any}}
+
+General container to sort the boundary conditions by type for the unstructured quadrilateral solver.
+It stores a set of global indices for each boundary condition type to expedite computation
+during the call to `calc_boundary_flux!`. The original dictionary form of the boundary conditions
+set by the user in the elixir file is also stored for printing.
+
+!!! warning "Experimental code"
+    This boundary condition container is experimental and can change any time.
+"""
+struct UnstructuredQuadSortedBoundaryTypes{N, BCs<:NTuple{N, Any}}
+  boundary_condition_types::BCs # specific boundary condition type(s), e.g. BoundaryConditionWall
+  boundary_indices::NTuple{N, Vector{Int}} # integer vectors containing global boundary indices
+  boundary_dictionary::Dict{Symbol, Any} # boundary conditions as set by the user in the elixir file
+end
+
+
+# constructor that "eats" the original boundary condition dictionary and sorts the information
+# from the `UnstructuredBoundaryContainer2D` in cache.boundaries according to the boundary types
+# and stores the associated global boundary indexing in NTuple
+function UnstructuredQuadSortedBoundaryTypes(boundary_conditions::Dict, cache)
+
+  # extract the unique boundary function routines from the dictionary
+  boundary_condition_types = Tuple(unique(collect(values(boundary_conditions))))
+  n_boundary_types = length(boundary_condition_types)
+
+  # pull and sort the indexing for each boundary type
+  _boundary_indices = Vector{Any}(nothing, n_boundary_types)
+  for j in 1:n_boundary_types
+    indices_for_current_type = Int[]
+    for (test_name, test_condition) in boundary_conditions
+      temp_indices = findall(x->x===test_name, cache.boundaries.name)
+      if test_condition === boundary_condition_types[j]
+        indices_for_current_type = vcat(indices_for_current_type, temp_indices)
+      end
+    end
+    _boundary_indices[j] = sort!(indices_for_current_type)
+  end
+
+  # convert the work array with the boundary indices into a tuple
+  boundary_indices = Tuple(_boundary_indices)
+
+  boundary_dictionary = boundary_conditions
+
+  return UnstructuredQuadSortedBoundaryTypes{n_boundary_types, typeof(boundary_condition_types)}(
+                                                                         boundary_condition_types,
+                                                                         boundary_indices,
+                                                                         boundary_dictionary)
+end