Add flexibility to nonconservative BCs

examples/dgmulti_2d/elixir_mhd_reflective_wall.jl

@@ -49,9 +49,9 @@ mesh = DGMultiMesh(solver, cells_per_dimension; periodicity = (false, false),
 # Create a "reflective-like" boundary condition by mirroring the velocity but leaving the magnetic field alone.
 # Note that this boundary condition is probably not entropy stable.
 function boundary_condition_velocity_slip_wall(u_inner, normal_direction::AbstractVector,
-                                               x, t, surface_flux_function,
+                                               x, t, surface_flux_functions,
                                                equations::IdealGlmMhdEquations2D)
-
+    surface_flux_function, nonconservative_flux_function = surface_flux_functions
     # Normalize the vector without using `normalize` since we need to multiply by the `norm_` later
     norm_ = norm(normal_direction)
     normal = normal_direction / norm_
@@ -63,8 +63,9 @@ function boundary_condition_velocity_slip_wall(u_inner, normal_direction::Abstra
     u_mirror = prim2cons(SVector(rho, v1 - 2 * v_normal * normal[1],
                                  v2 - 2 * v_normal * normal[2],
                                  v3, p, B1, B2, B3, psi), equations)
-
-    return surface_flux_function(u_inner, u_mirror, normal, equations) * norm_
+    flux = surface_flux_function(u_inner, u_mirror, normal, equations) * norm_
+    noncons_flux = nonconservative_flux_function(u_inner, u_mirror, normal, equations) * norm_
+    return flux, noncons_flux
 end
 
 boundary_conditions = (; x_neg = boundary_condition_velocity_slip_wall,

src/basic_types.jl

@@ -74,13 +74,37 @@ struct BoundaryConditionDoNothing end
     return surface_flux(u_inner, u_inner, orientation_or_normal_direction, equations)
 end
 
+# This version can be called by hyperbolic solvers on logically Cartesian meshes
+@inline function (::BoundaryConditionDoNothing)(u_inner,
+                                                orientation_or_normal_direction,
+                                                direction::Integer, x, t,
+                                                surface_flux_functions::Tuple,
+                                                equations)
+    surface_flux_function, nonconservative_flux_function = surface_flux_functions
+    return surface_flux_function(u_inner, u_inner, orientation_or_normal_direction,
+                                 equations),
+           nonconservative_flux_function(u_inner, u_inner,
+                                         orientation_or_normal_direction, equations)
+end
+
 # This version can be called by hyperbolic solvers on unstructured, curved meshes
 @inline function (::BoundaryConditionDoNothing)(u_inner,
                                                 outward_direction::AbstractVector,
                                                 x, t, surface_flux, equations)
     return surface_flux(u_inner, u_inner, outward_direction, equations)
 end
 
+@inline function (::BoundaryConditionDoNothing)(u_inner,
+                                                outward_direction::AbstractVector,
+                                                x, t, surface_flux_functions::Tuple,
+                                                equations)
+    surface_flux_function, nonconservative_flux_function = surface_flux_functions
+
+    return surface_flux_function(u_inner, u_inner, outward_direction, equations),
+           nonconservative_flux_function(u_inner, u_inner, outward_direction,
+                                         equations)
+end
+
 # This version can be called by parabolic solvers
 @inline function (::BoundaryConditionDoNothing)(inner_flux_or_state, other_args...)
     return inner_flux_or_state

src/equations/equations.jl

@@ -181,6 +181,35 @@ end
     return flux
 end
 
+# Dirichlet-type boundary condition for use with TreeMesh or StructuredMesh
+@inline function (boundary_condition::BoundaryConditionDirichlet)(u_inner,
+                                                                  orientation_or_normal,
+                                                                  direction,
+                                                                  x, t,
+                                                                  surface_flux_functions::Tuple,
+                                                                  equations)
+    surface_flux_function, nonconservative_flux_function = surface_flux_functions
+
+    u_boundary = boundary_condition.boundary_value_function(x, t, equations)
+
+    # Calculate boundary flux
+    if iseven(direction) # u_inner is "left" of boundary, u_boundary is "right" of boundary
+        flux = surface_flux_function(u_inner, u_boundary, orientation_or_normal,
+                                     equations)
+        noncons_flux = nonconservative_flux_function(u_inner, u_boundary,
+                                                orientation_or_normal,
+                                                equations)
+    else # u_boundary is "left" of boundary, u_inner is "right" of boundary
+        flux = surface_flux_function(u_boundary, u_inner, orientation_or_normal,
+                                     equations)
+        noncons_flux = nonconservative_flux_function(u_boundary, u_inner,
+                                                orientation_or_normal,
+                                                equations)
+    end
+
+    return flux, noncons_flux
+end
+
 # Dirichlet-type boundary condition for use with UnstructuredMesh2D
 # Note: For unstructured we lose the concept of an "absolute direction"
 @inline function (boundary_condition::BoundaryConditionDirichlet)(u_inner,
@@ -197,6 +226,25 @@ end
     return flux
 end
 
+# Dirichlet-type boundary condition for equations with non-conservative terms for use with UnstructuredMesh2D
+# Note: For unstructured we lose the concept of an "absolute direction"
+@inline function (boundary_condition::BoundaryConditionDirichlet)(u_inner,
+                                                                  normal_direction::AbstractVector,
+                                                                  x, t,
+                                                                  surface_flux_functions::Tuple,
+                                                                  equations)
+    surface_flux_function, nonconservative_flux_function = surface_flux_functions
+
+    # get the external value of the solution
+    u_boundary = boundary_condition.boundary_value_function(x, t, equations)
+
+    # Calculate boundary flux
+    flux = surface_flux_function(u_inner, u_boundary, normal_direction, equations)
+    noncons_flux = nonconservative_flux_function(u_inner, u_boundary, normal_direction,
+                                            equations)
+    return flux, noncons_flux
+end
+
 # operator types used for dispatch on parabolic boundary fluxes
 struct Gradient end
 struct Divergence end

src/equations/shallow_water_1d.jl

@@ -160,9 +160,11 @@ For details see Section 9.2.5 of the book:
 """
 @inline function boundary_condition_slip_wall(u_inner, orientation_or_normal, direction,
                                               x, t,
-                                              surface_flux_function,
+                                              surface_flux_functions,
                                               equations::ShallowWaterEquations1D)
 
+    # The boundary conditions for the non-conservative term are identically 0 here.
+    surface_flux_function, nonconservative_flux_function = surface_flux_functions
     # create the "external" boundary solution state
     u_boundary = SVector(u_inner[1],
                          -u_inner[2],
@@ -172,12 +174,18 @@ For details see Section 9.2.5 of the book:
     if iseven(direction) # u_inner is "left" of boundary, u_boundary is "right" of boundary
         flux = surface_flux_function(u_inner, u_boundary, orientation_or_normal,
                                      equations)
+        noncons_flux = nonconservative_flux_function(u_inner, u_boundary,
+                                                orientation_or_normal,
+                                                equations)
     else # u_boundary is "left" of boundary, u_inner is "right" of boundary
         flux = surface_flux_function(u_boundary, u_inner, orientation_or_normal,
                                      equations)
+        noncons_flux = nonconservative_flux_function(u_boundary, u_inner,
+                                                orientation_or_normal,
+                                                equations)
     end
 
-    return flux
+    return flux, noncons_flux
 end
 
 # Calculate 1D flux for a single point

src/equations/shallow_water_2d.jl

@@ -180,8 +180,10 @@ For details see Section 9.2.5 of the book:
 """
 @inline function boundary_condition_slip_wall(u_inner, normal_direction::AbstractVector,
                                               x, t,
-                                              surface_flux_function,
+                                              surface_flux_functions,
                                               equations::ShallowWaterEquations2D)
+    surface_flux_function, nonconservative_flux_function = surface_flux_functions
+
     # normalize the outward pointing direction
     normal = normal_direction / norm(normal_direction)
 
@@ -196,8 +198,10 @@ For details see Section 9.2.5 of the book:
 
     # calculate the boundary flux
     flux = surface_flux_function(u_inner, u_boundary, normal_direction, equations)
+    noncons_flux = nonconservative_flux_function(u_inner, u_boundary, normal_direction,
+                                            equations)
 
-    return flux
+    return flux, noncons_flux
 end
 
 """
@@ -208,8 +212,10 @@ Should be used together with [`TreeMesh`](@ref).
 """
 @inline function boundary_condition_slip_wall(u_inner, orientation,
                                               direction, x, t,
-                                              surface_flux_function,
+                                              surface_flux_functions,
                                               equations::ShallowWaterEquations2D)
+    # The boundary conditions for the non-conservative term are identically 0 here.
+    surface_flux_function, nonconservative_flux_function = surface_flux_functions
     ## get the appropriate normal vector from the orientation
     if orientation == 1
         u_boundary = SVector(u_inner[1], -u_inner[2], u_inner[3], u_inner[4])
@@ -220,11 +226,15 @@ Should be used together with [`TreeMesh`](@ref).
     # Calculate boundary flux
     if iseven(direction) # u_inner is "left" of boundary, u_boundary is "right" of boundary
         flux = surface_flux_function(u_inner, u_boundary, orientation, equations)
+        noncons_flux = nonconservative_flux_function(u_inner, u_boundary, orientation,
+                                                equations)
     else # u_boundary is "left" of boundary, u_inner is "right" of boundary
         flux = surface_flux_function(u_boundary, u_inner, orientation, equations)
+        noncons_flux = nonconservative_flux_function(u_boundary, u_inner, orientation,
+                                                equations)
     end
 
-    return flux
+    return flux, noncons_flux
 end
 
 # Calculate 1D flux for a single point

src/solvers/dgmulti/dg.jl

@@ -519,7 +519,6 @@ function calc_single_boundary_flux!(cache, t, boundary_condition, boundary_key,
                                     dg::DGMulti{NDIMS}) where {NDIMS}
     rd = dg.basis
     md = mesh.md
-    surface_flux, nonconservative_flux = dg.surface_integral.surface_flux
 
     # reshape face/normal arrays to have size = (num_points_on_face, num_faces_total).
     # mesh.boundary_faces indexes into the columns of these face-reshaped arrays.
@@ -546,21 +545,16 @@ function calc_single_boundary_flux!(cache, t, boundary_condition, boundary_key,
 
             # Compute conservative and non-conservative fluxes separately.
             # This imposes boundary conditions on the conservative part of the flux.
-            cons_flux_at_face_node = boundary_condition(u_face_values[i, f],
-                                                        face_normal, face_coordinates,
-                                                        t,
-                                                        surface_flux, equations)
-
-            # Compute pointwise nonconservative numerical flux at the boundary.
-            noncons_flux_at_face_node = boundary_condition(u_face_values[i, f],
-                                                           face_normal,
-                                                           face_coordinates,
-                                                           t,
-                                                           nonconservative_flux,
-                                                           equations)
+            cons_flux_at_face_node, noncons_flux_at_face_node = boundary_condition(u_face_values[i,
+                                                                                                 f],
+                                                                                   face_normal,
+                                                                                   face_coordinates,
+                                                                                   t,
+                                                                                   dg.surface_integral.surface_flux,
+                                                                                   equations)
 
             flux_face_values[i, f] = (cons_flux_at_face_node +
-                                      0.5 * noncons_flux_at_face_node) * Jf[i, f]
+                                      0.5f0 * noncons_flux_at_face_node) * Jf[i, f]
         end
     end
 

src/solvers/dgsem_p4est/dg_2d.jl

@@ -349,7 +349,6 @@ end
                                      boundary_index)
     @unpack boundaries = cache
     @unpack node_coordinates, contravariant_vectors = cache.elements
-    surface_flux, nonconservative_flux = surface_integral.surface_flux
 
     # Extract solution data from boundary container
     u_inner = get_node_vars(boundaries.u, equations, dg, node_index, boundary_index)
@@ -364,11 +363,8 @@ end
 
     # Call pointwise numerical flux function for the conservative part
     # in the normal direction on the boundary
-    flux_ = boundary_condition(u_inner, normal_direction, x, t, surface_flux, equations)
-
-    # Compute pointwise nonconservative numerical flux at the boundary.
-    noncons_ = boundary_condition(u_inner, normal_direction, x, t, nonconservative_flux,
-                                  equations)
+    flux, noncons_flux = boundary_condition(u_inner, normal_direction, x, t,
+                                         surface_integral.surface_flux, equations)
 
     # Copy flux to element storage in the correct orientation
     for v in eachvariable(equations)

src/solvers/dgsem_tree/dg_1d.jl

@@ -564,7 +564,6 @@ function calc_boundary_flux_by_direction!(surface_flux_values::AbstractArray{<:A
                                           nonconservative_terms::True, equations,
                                           surface_integral, dg::DG, cache,
                                           direction, first_boundary, last_boundary)
-    surface_flux, nonconservative_flux = surface_integral.surface_flux
     @unpack u, neighbor_ids, neighbor_sides, node_coordinates, orientations = cache.boundaries
 
     @threaded for boundary in first_boundary:last_boundary
@@ -579,12 +578,11 @@ function calc_boundary_flux_by_direction!(surface_flux_values::AbstractArray{<:A
             u_inner = u_rr
         end
         x = get_node_coords(node_coordinates, equations, dg, boundary)
-        flux = boundary_condition(u_inner, orientations[boundary], direction, x, t,
-                                  surface_flux,
-                                  equations)
-        noncons_flux = boundary_condition(u_inner, orientations[boundary], direction, x,
-                                          t, nonconservative_flux,
-                                          equations)
+
+        flux, noncons_flux = boundary_condition(u_inner, orientations[boundary],
+                                                direction, x, t,
+                                                surface_integral.surface_flux,
+                                                equations)
 
         # Copy flux to left and right element storage
         for v in eachvariable(equations)

src/solvers/dgsem_tree/dg_2d.jl

@@ -762,7 +762,6 @@ function calc_boundary_flux_by_direction!(surface_flux_values::AbstractArray{<:A
                                           nonconservative_terms::True, equations,
                                           surface_integral, dg::DG, cache,
                                           direction, first_boundary, last_boundary)
-    surface_flux, nonconservative_flux = surface_integral.surface_flux
     @unpack u, neighbor_ids, neighbor_sides, node_coordinates, orientations = cache.boundaries
 
     @threaded for boundary in first_boundary:last_boundary
@@ -778,12 +777,10 @@ function calc_boundary_flux_by_direction!(surface_flux_values::AbstractArray{<:A
                 u_inner = u_rr
             end
             x = get_node_coords(node_coordinates, equations, dg, i, boundary)
-            flux = boundary_condition(u_inner, orientations[boundary], direction, x, t,
-                                      surface_flux,
-                                      equations)
-            noncons_flux = boundary_condition(u_inner, orientations[boundary],
-                                              direction, x, t, nonconservative_flux,
-                                              equations)
+            flux, noncons_flux = boundary_condition(u_inner, orientations[boundary],
+                                                    direction, x, t,
+                                                    surface_integral.surface_flux,
+                                                    equations)
 
             # Copy flux to left and right element storage
             for v in eachvariable(equations)

src/solvers/dgsem_unstructured/dg_2d.jl

@@ -426,7 +426,6 @@ end
                                      surface_integral, dg::DG, cache,
                                      node_index, side_index, element_index,
                                      boundary_index)
-    surface_flux, nonconservative_flux = surface_integral.surface_flux
     @unpack normal_directions = cache.elements
     @unpack u, node_coordinates = cache.boundaries
 
@@ -442,11 +441,10 @@ end
 
     # Call pointwise numerical flux function for the conservative part
     # in the normal direction on the boundary
-    flux = boundary_condition(u_inner, outward_direction, x, t, surface_flux, equations)
+    flux, noncons_flux = boundary_condition(u_inner, outward_direction, x, t,
+                                            surface_integral.surface_flux, equations)
 
     # Compute pointwise nonconservative numerical flux at the boundary.
-    noncons_flux = boundary_condition(u_inner, outward_direction, x, t,
-                                      nonconservative_flux, equations)
 
     for v in eachvariable(equations)
         # Note the factor 0.5 necessary for the nonconservative fluxes based on

test/test_type.jl

@@ -2258,7 +2258,7 @@ isdir(outdir) && rm(outdir, recursive = true)
             directions = [1, 2]
             normal_direction = SVector(one(RealT))
 
-            surface_flux_function = flux_lax_friedrichs
+            surface_flux_function = (flux_lax_friedrichs, flux_lax_friedrichs)
             dissipation = DissipationLocalLaxFriedrichs()
             numflux = FluxHLL()
 
@@ -2347,7 +2347,7 @@ isdir(outdir) && rm(outdir, recursive = true)
             directions = [1, 2, 3, 4]
             normal_direction = SVector(one(RealT), zero(RealT))
 
-            surface_flux_function = flux_lax_friedrichs
+            surface_flux_function = (flux_lax_friedrichs, flux_lax_friedrichs)
             dissipation = DissipationLocalLaxFriedrichs()
             numflux = FluxHLL()