unstructured BCs (#573)

* improve uniform_flow_state * use symbols instead of strings for boundary names * introduce a function barrier for boundary fluxes * fix unstructured periodic BCs Co-authored-by: Michael Schlottke-Lakemper <[email protected]>
trixi-framework · May 12, 2021 · ca0a78f · ca0a78f
1 parent 0d1bf98
commit ca0a78f
Show file tree

Hide file tree

Showing 6 changed files with 45 additions and 37 deletions.
diff --git a/examples/2d/elixir_euler_unstructured_quad_basic.jl b/examples/2d/elixir_euler_unstructured_quad_basic.jl
@@ -12,11 +12,11 @@ initial_condition = initial_condition_convergence_test
 source_terms = source_terms_convergence_test
 
 boundary_condition_convergence_test = BoundaryConditionDirichlet(initial_condition)
-boundary_conditions = Dict( "Slant" => boundary_condition_convergence_test,
-                            "Bezier" => boundary_condition_convergence_test,
-                            "Right"  => boundary_condition_convergence_test,
-                            "Bottom" => boundary_condition_convergence_test,
-                            "Top"    => boundary_condition_convergence_test )
+boundary_conditions = Dict( :Slant  => boundary_condition_convergence_test,
+                            :Bezier => boundary_condition_convergence_test,
+                            :Right  => boundary_condition_convergence_test,
+                            :Bottom => boundary_condition_convergence_test,
+                            :Top    => boundary_condition_convergence_test )
 
 ###############################################################################
 # Get the DG approximation space

diff --git a/examples/2d/elixir_euler_unstructured_quad_free_stream.jl b/examples/2d/elixir_euler_unstructured_quad_free_stream.jl
@@ -11,13 +11,13 @@ equations = CompressibleEulerEquations2D(1.4)
 initial_condition = initial_condition_constant
 
 boundary_condition_free_stream = BoundaryConditionDirichlet(initial_condition)
-boundary_conditions = Dict( "Body"    => boundary_condition_free_stream,
-                            "Button1" => boundary_condition_free_stream,
-                            "Button2" => boundary_condition_free_stream,
-                            "Eye1"    => boundary_condition_free_stream,
-                            "Eye2"    => boundary_condition_free_stream,
-                            "Smile"   => boundary_condition_free_stream,
-                            "Bowtie"  => boundary_condition_free_stream )
+boundary_conditions = Dict( :Body    => boundary_condition_free_stream,
+                            :Button1 => boundary_condition_free_stream,
+                            :Button2 => boundary_condition_free_stream,
+                            :Eye1    => boundary_condition_free_stream,
+                            :Eye2    => boundary_condition_free_stream,
+                            :Smile   => boundary_condition_free_stream,
+                            :Bowtie  => boundary_condition_free_stream )
 
 ###############################################################################
 # Get the DG approximation space

diff --git a/examples/2d/elixir_euler_unstructured_quad_wall_bc.jl b/examples/2d/elixir_euler_unstructured_quad_wall_bc.jl
@@ -8,16 +8,17 @@ using Trixi
 
 equations = CompressibleEulerEquations2D(1.4)
 
-function uniform_flow_state(x, t, equations)
+@inline function uniform_flow_state(x, t, equations::CompressibleEulerEquations2D)
 
   # set the freestream flow parameters
   rho_freestream = 1.0
   u_freestream = 0.3
   p_freestream = inv(equations.gamma)
 
   theta = pi / 90.0 # analogous with a two degree angle of attack
-  v1 = u_freestream * cos(theta)
-  v2 = u_freestream * sin(theta)
+  si, co = sincos(theta)
+  v1 = u_freestream * co
+  v2 = u_freestream * si
 
   prim = SVector(rho_freestream, v1, v2, p_freestream)
   return prim2cons(prim, equations)
@@ -27,11 +28,11 @@ initial_condition = uniform_flow_state
 
 boundary_condition_uniform_flow = BoundaryConditionDirichlet(uniform_flow_state)
 boundary_condition_slip_wall = BoundaryConditionWall(boundary_state_slip_wall)
-boundary_conditions = Dict( "Bottom" => boundary_condition_uniform_flow,
-                            "Top"    => boundary_condition_uniform_flow,
-                            "Right"  => boundary_condition_uniform_flow,
-                            "Left"   => boundary_condition_uniform_flow,
-                            "Circle" => boundary_condition_slip_wall )
+boundary_conditions = Dict( :Bottom => boundary_condition_uniform_flow,
+                            :Top    => boundary_condition_uniform_flow,
+                            :Right  => boundary_condition_uniform_flow,
+                            :Left   => boundary_condition_uniform_flow,
+                            :Circle => boundary_condition_slip_wall )
 
 ###############################################################################
 # Get the DG approximation space

diff --git a/src/mesh/unstructured_quad_mesh.jl b/src/mesh/unstructured_quad_mesh.jl
@@ -21,7 +21,7 @@ mutable struct UnstructuredQuadMesh{RealT<:Real, CurvedSurfaceT<:CurvedSurface{R
   polydeg              ::Int
   corners              ::Array{RealT, 2}  # [ndims, n_corners]
   neighbour_information::Array{Int, 2}  # [neighbour node/element/edge ids, n_surfaces]
-  boundary_names       ::Array{String, 2} # [local sides, n_elements]
+  boundary_names       ::Array{Symbol, 2} # [local sides, n_elements]
   periodicity          ::Bool
   element_node_ids     ::Array{Int, 2} # [node ids, n_elements]
   element_is_curved    ::Vector{Bool}
@@ -62,7 +62,7 @@ function UnstructuredQuadMesh(filename; RealT=Float64, periodicity=false, unsave
   element_is_curved = Array{Bool}(undef, n_elements)
   CurvedSurfaceT    = CurvedSurface{RealT}
   surface_curves    = Array{CurvedSurfaceT}(undef, (4, n_elements))
-  bndy_names        = Array{String}(undef, (4, n_elements))
+  boundary_names    = Array{Symbol}(undef, (4, n_elements))
 
   # create the Chebyshev-Gauss-Lobatto nodes used to represent any curved boundaries that are
   # required to construct the sides
@@ -73,7 +73,7 @@ function UnstructuredQuadMesh(filename; RealT=Float64, periodicity=false, unsave
 
   arrays = (; corner_nodes, interface_info, element_node_ids, curved_check,
               cornerNodeVals, tempNodes, curve_vals,
-              element_is_curved, surface_curves, bndy_names)
+              element_is_curved, surface_curves, boundary_names)
   counters = (; n_corners, n_surfaces, n_elements)
 
   n_boundaries = parse_mesh_file!(arrays, RealT, CurvedSurfaceT, file_lines, counters, cheby_nodes, bary_weights)
@@ -89,14 +89,14 @@ function UnstructuredQuadMesh(filename; RealT=Float64, periodicity=false, unsave
   return UnstructuredQuadMesh{RealT, CurvedSurfaceT}(
     filename, n_corners, n_surfaces, n_interfaces, n_boundaries,
     n_elements, mesh_polydeg, corner_nodes,
-    interface_info, bndy_names, periodicity,
+    interface_info, boundary_names, periodicity,
     element_node_ids, element_is_curved, surface_curves, "", unsaved_changes)
 end
 
 function parse_mesh_file!(arrays, RealT, CurvedSurfaceT, file_lines, counters, cheby_nodes, bary_weights)
   @unpack ( corner_nodes, interface_info, element_node_ids, curved_check,
             cornerNodeVals, tempNodes, curve_vals,
-            element_is_curved, surface_curves, bndy_names ) = arrays
+            element_is_curved, surface_curves, boundary_names ) = arrays
   @unpack n_corners, n_surfaces, n_elements = counters
   mesh_nnodes = length(cheby_nodes)
 
@@ -166,7 +166,7 @@ function parse_mesh_file!(arrays, RealT, CurvedSurfaceT, file_lines, counters, c
       element_is_curved[j] = false
       file_idx  += 1
       # read all the boundary names
-      bndy_names[:, j] = split(file_lines[file_idx])
+      boundary_names[:, j] = map(Symbol, split(file_lines[file_idx]))
     else
       # quadrilateral element has at least one curved side
       element_is_curved[j] = true
@@ -211,7 +211,7 @@ function parse_mesh_file!(arrays, RealT, CurvedSurfaceT, file_lines, counters, c
       end
       # finally read in the boundary names where "---" means an internal connection
       file_idx  += 1
-      bndy_names[:, j] = split(file_lines[file_idx])
+      boundary_names[:, j] = map(Symbol, split(file_lines[file_idx]))
     end
     # one last increment to the global index to read the next piece of element information
     file_idx += 1

diff --git a/src/solvers/dg_unstructured_quad/containers_2d.jl b/src/solvers/dg_unstructured_quad/containers_2d.jl
@@ -195,16 +195,16 @@ function init_interfaces!(interfaces, edge_information, boundary_names, polydeg,
       primary_element = edge_information[3,j]
       # Note: This is a way to get the neighbour element number and local side from a square
       #       structured mesh where the element local surface numbering is right-handed
-      if boundary_names[primary_side, primary_element] == "Bottom"
+      if boundary_names[primary_side, primary_element] === :Bottom
         secondary_element = primary_element + (n_elements - convert(Int, sqrt(n_elements)))
         secondary_side    = 3
-      elseif boundary_names[primary_side, primary_element] == "Top"
+      elseif boundary_names[primary_side, primary_element] === :Top
         secondary_element = primary_element - (n_elements - convert(Int, sqrt(n_elements)))
         secondary_side    = 1
-      elseif boundary_names[primary_side, primary_element] == "Left"
+      elseif boundary_names[primary_side, primary_element] === :Left
         secondary_element = primary_element + (convert(Int, sqrt(n_elements)) - 1)
         secondary_side    = 2
-      elseif boundary_names[primary_side, primary_element] == "Right"
+      elseif boundary_names[primary_side, primary_element] === :Right
         secondary_element = primary_element - (convert(Int, sqrt(n_elements)) - 1)
         secondary_side    = 4
       end
@@ -229,7 +229,7 @@ struct UnstructuredBoundaryContainer2D{RealT<:Real, uEltype<:Real, NVARS, POLYDE
   element_id      ::Vector{Int}       # [boundaries]
   element_side_id ::Vector{Int}       # [boundaries]
   node_coordinates::Array{RealT, 3}   # [ndims, nnodes, boundaries]
-  name            ::Vector{String}    # [boundaries]
+  name            ::Vector{Symbol}    # [boundaries]
 end
 
 
@@ -245,7 +245,7 @@ function UnstructuredBoundaryContainer2D{RealT, uEltype, NVARS, POLYDEG}(capacit
   element_id       = fill(typemin(Int), capacity)
   element_side_id  = fill(typemin(Int), capacity)
   node_coordinates = fill(nan_RealT, (2, n_nodes, capacity))
-  name             = fill("empty", capacity)
+  name             = fill(:empty, capacity)
 
   return UnstructuredBoundaryContainer2D{RealT, uEltype, NVARS, POLYDEG}(u, element_id, element_side_id,
                                                                          node_coordinates, name)

diff --git a/src/solvers/dg_unstructured_quad/dg_2d.jl b/src/solvers/dg_unstructured_quad/dg_2d.jl
@@ -235,15 +235,22 @@ function calc_boundary_flux!(cache, t, boundary_conditions,
     boundary_condition = boundary_conditions[ name[boundary] ]
 
     # calc boundary flux on the current boundary interface
-    for j in eachnode(dg)
-      calc_boundary_flux!(surface_flux_values, t, boundary_condition, mesh, equations, dg, cache,
-                          j, side, element, boundary)
-    end
+    calc_boundary_flux!(surface_flux_values, t, boundary_condition, mesh, equations, dg, cache,
+                        side, element, boundary)
   end
 
   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)
+  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