WIP: Non-Conservative BCs update

examples/unstructured_2d_dgsem/elixir_shallowwater_three_mound_dam_break.jl

@@ -53,15 +53,17 @@
 initial_condition = initial_condition_three_mounds
 
 function boundary_condition_outflow(u_inner, normal_direction::AbstractVector, x, t,
-                                    surface_flux_function,
+                                    surface_flux_functions,
                                     equations::ShallowWaterEquationsWetDry2D)
+    surface_flux_function, nonconservative_flux_function = surface_flux_functions
     # Impulse and bottom from inside, height from external state
     u_outer = SVector(equations.threshold_wet, u_inner[2], u_inner[3], u_inner[4])
 
     # calculate the boundary flux
     flux = surface_flux_function(u_inner, u_outer, normal_direction, equations)
-
-    return flux
+    noncons_flux = nonconservative_flux_function(u_inner, u_outer, normal_direction,
+                                                 equations)
+    return flux, noncons_flux
 end
 
 boundary_conditions = Dict(:Bottom => boundary_condition_slip_wall,

src/equations/shallow_water_exner_1d.jl

@@ -190,8 +190,10 @@
 @inline function Trixi.boundary_condition_slip_wall(u_inner, orientation_or_normal,
                                                     direction,
                                                     x, t,
-                                                    surface_flux_function,
+                                                    surface_flux_functions,
                                                     equations::ShallowWaterExnerEquations1D)
+    surface_flux_function, nonconservative_flux_function = surface_flux_functions
+
     # create the "external" boundary solution state
     u_boundary = SVector(u_inner[1],
                          -u_inner[2],
@@ -201,12 +203,18 @@
     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
 
 # Set initial conditions at physical location `x` for time `t`

src/equations/shallow_water_multilayer_1d.jl

@@ -256,8 +256,10 @@
 """
 @inline function Trixi.boundary_condition_slip_wall(u_inner, orientation_or_normal,
                                                     direction,
-                                                    x, t, surface_flux_function,
+                                                    x, t, surface_flux_functions,
                                                     equations::ShallowWaterMultiLayerEquations1D)
+    surface_flux_function, nonconservative_flux_function = surface_flux_functions
+
     # Create the "external" boundary solution state
     h = waterheight(u_inner, equations)
     hv = momentum(u_inner, equations)
@@ -267,11 +269,19 @@
 
     # Calculate the boundary flux
     if iseven(direction) # u_inner is "left" of boundary, u_boundary is "right" of boundary
-        f = surface_flux_function(u_inner, u_boundary, orientation_or_normal, equations)
+        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
-        f = surface_flux_function(u_boundary, u_inner, orientation_or_normal, equations)
+        flux = surface_flux_function(u_boundary, u_inner, orientation_or_normal,
+                                     equations)
+        noncons_flux = nonconservative_flux_function(u_inner, u_boundary,
+                                                     orientation_or_normal,
+                                                     equations)
     end
-    return f
+    return flux, noncons_flux
 end
 
 # Calculate 1D advective portion of the flux for a single point

src/equations/shallow_water_multilayer_2d.jl

@@ -331,8 +331,10 @@
 @inline function Trixi.boundary_condition_slip_wall(u_inner,
                                                     normal_direction::AbstractVector,
                                                     x, t,
-                                                    surface_flux_function,
+                                                    surface_flux_functions,
                                                     equations::ShallowWaterMultiLayerEquations2D)
+    surface_flux_function, nonconservative_flux_function = surface_flux_functions
+
     # normalize the outward pointing direction
     normal = normal_direction / norm(normal_direction)
 
@@ -352,8 +354,10 @@
 
     # 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
 
 """
@@ -364,8 +368,10 @@
 """
 @inline function Trixi.boundary_condition_slip_wall(u_inner, orientation,
                                                     direction, x, t,
-                                                    surface_flux_function,
+                                                    surface_flux_functions,
                                                     equations::ShallowWaterMultiLayerEquations2D)
+    surface_flux_function, nonconservative_flux_function = surface_flux_functions
+
     # Extract internal values
     h = waterheight(u_inner, equations)
     h_v1, h_v2 = momentum(u_inner, equations)
@@ -381,11 +387,15 @@
     # 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 2D advective portion of the flux for a single point

src/equations/shallow_water_two_layer_1d.jl

@@ -174,8 +174,10 @@
 """
 @inline function Trixi.boundary_condition_slip_wall(u_inner, orientation_or_normal,
                                                     direction,
-                                                    x, t, surface_flux_function,
+                                                    x, t, surface_flux_functions,
                                                     equations::ShallowWaterTwoLayerEquations1D)
+    surface_flux_function, nonconservative_flux_function = surface_flux_functions
+
     # create the "external" boundary solution state
     u_boundary = SVector(u_inner[1],
                          -u_inner[2],
@@ -185,11 +187,18 @@
 
     # calculate the boundary flux
     if iseven(direction) # u_inner is "left" of boundary, u_boundary is "right" of boundary
-        f = surface_flux_function(u_inner, u_boundary, orientation_or_normal, equations)
+        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
-        f = surface_flux_function(u_boundary, u_inner, orientation_or_normal, equations)
+        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 f
+
+    return flux, noncons_flux
 end
 
 # Calculate 1D flux for a single point

src/equations/shallow_water_two_layer_2d.jl

@@ -238,8 +238,10 @@
 """
 @inline function Trixi.boundary_condition_slip_wall(u_inner,
                                                     normal_direction::AbstractVector,
-                                                    x, t, surface_flux_function,
+                                                    x, t, surface_flux_functions,
                                                     equations::ShallowWaterTwoLayerEquations2D)
+    surface_flux_function, nonconservative_flux_function = surface_flux_functions
+
     # normalize the outward pointing direction
     normal = normal_direction / norm(normal_direction)
 
@@ -258,7 +260,10 @@
 
     # calculate the boundary flux
     flux = surface_flux_function(u_inner, u_boundary, normal_direction, equations)
-    return flux
+    noncons_flux = nonconservative_flux_function(u_inner, u_boundary, normal_direction,
+                                                 equations)
+
+    return flux, noncons_flux
 end
 
 # Calculate 1D flux for a single point

src/equations/shallow_water_wet_dry_1d.jl

@@ -179,8 +179,10 @@
 @inline function Trixi.boundary_condition_slip_wall(u_inner, orientation_or_normal,
                                                     direction,
                                                     x, t,
-                                                    surface_flux_function,
+                                                    surface_flux_functions,
                                                     equations::ShallowWaterEquationsWetDry1D)
+    surface_flux_function, nonconservative_flux_function = surface_flux_functions
+
     # This can not be dispatched, when Flux Hydrostactic reconstruction is used
     # create the "external" boundary solution state
     u_boundary = SVector(u_inner[1],
@@ -191,12 +193,18 @@
     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_wet_dry_2d.jl

@@ -177,8 +177,10 @@
 @inline function Trixi.boundary_condition_slip_wall(u_inner,
                                                     normal_direction::AbstractVector,
                                                     x, t,
-                                                    surface_flux_function,
+                                                    surface_flux_functions,
                                                     equations::ShallowWaterEquationsWetDry2D)
+    surface_flux_function, nonconservative_flux_function = surface_flux_functions
+
     # normalize the outward pointing direction
     normal = normal_direction / norm(normal_direction)
 
@@ -193,8 +195,9 @@
 
     # calculate the boundary flux
     flux = surface_flux_function(u_inner, u_boundary, normal_direction, equations)
-
-    return flux
+    noncons_flux = nonconservative_flux_function(u_inner, u_boundary, normal_direction,
+                                                 equations)
+    return flux, noncons_flux
 end
 
 """
@@ -205,8 +208,10 @@
 """
 @inline function Trixi.boundary_condition_slip_wall(u_inner, orientation,
                                                     direction, x, t,
-                                                    surface_flux_function,
+                                                    surface_flux_functions,
                                                     equations::ShallowWaterEquationsWetDry2D)
+    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])
@@ -217,11 +222,15 @@
     # 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