Remove normal_direction_ll for nonconservative terms

NEWS.md

@@ -6,13 +6,27 @@ for human readability.
 
 ## Changes when updating to v0.9 from v0.8.x
 
+#### Added
+
+- Boundary conditions are now supported on nonconservative terms ([#2062]).
+
 #### Changed
 
 - We removed the first argument `semi` corresponding to a `Semidiscretization` from the 
   `AnalysisSurfaceIntegral` constructor, as it is no longer needed (see [#1959]).
-  The `AnalysisSurfaceIntegral` now only takes the arguments `boundary_symbols` and `variable`. ([#2069])
+  The `AnalysisSurfaceIntegral` now only takes the arguments `boundary_symbols` and `variable`.
+  ([#2069])
 - In functions `rhs!`, `rhs_parabolic!`  we removed the unused argument `initial_condition`. ([#2037])
   Users should not be affected by this.
+- Nonconservative terms depend only on `normal_direction_average` instead of both 
+  `normal_direction_average` and `normal_direction_ll`, such that the function signature is now 
+  identical with conservative fluxes. This required a change of the `normal_direction` in
+  `flux_nonconservative_powell` ([#2062]).
+
+#### Deprecated
+
+#### Removed
+
 
 ## Changes in the v0.8 lifecycle
 

examples/structured_2d_dgsem/elixir_mhd_coupled.jl

@@ -31,14 +31,6 @@ equations = IdealGlmMhdEquations2D(gamma)
 
 cells_per_dimension = (32, 64)
 
-# Extend the definition of the non-conservative Powell flux functions.
-import Trixi.flux_nonconservative_powell
-function flux_nonconservative_powell(u_ll, u_rr,
-                                     normal_direction_ll::AbstractVector,
-                                     equations::IdealGlmMhdEquations2D)
-    flux_nonconservative_powell(u_ll, u_rr, normal_direction_ll, normal_direction_ll,
-                                equations)
-end
 volume_flux = (flux_hindenlang_gassner, flux_nonconservative_powell)
 solver = DGSEM(polydeg = 3,
                surface_flux = (flux_lax_friedrichs, flux_nonconservative_powell),

src/basic_types.jl

@@ -71,14 +71,14 @@ struct BoundaryConditionDoNothing end
                                                 orientation_or_normal_direction,
                                                 direction::Integer, x, t, surface_flux,
                                                 equations)
-    return flux(u_inner, orientation_or_normal_direction, equations)
+    return surface_flux(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 flux(u_inner, outward_direction, equations)
+    return surface_flux(u_inner, u_inner, outward_direction, equations)
 end
 
 # This version can be called by parabolic solvers

src/equations/ideal_glm_mhd_2d.jl

@@ -196,18 +196,15 @@ end
     flux_nonconservative_powell(u_ll, u_rr, orientation::Integer,
                                 equations::IdealGlmMhdEquations2D)
     flux_nonconservative_powell(u_ll, u_rr,
-                                normal_direction_ll     ::AbstractVector,
-                                normal_direction_average::AbstractVector,
+                                normal_direction::AbstractVector,
                                 equations::IdealGlmMhdEquations2D)
 
 Non-symmetric two-point flux discretizing the nonconservative (source) term of
 Powell and the Galilean nonconservative term associated with the GLM multiplier
 of the [`IdealGlmMhdEquations2D`](@ref).
 
-On curvilinear meshes, this nonconservative flux depends on both the
-contravariant vector (normal direction) at the current node and the averaged
-one. This is different from numerical fluxes used to discretize conservative
-terms.
+On curvilinear meshes, the implementation differs from the reference since we use the averaged 
+normal direction for the metrics dealiasing.
 
 ## References
 - Marvin Bohm, Andrew R.Winters, Gregor J. Gassner, Dominik Derigs,
@@ -254,8 +251,7 @@ terms.
 end
 
 @inline function flux_nonconservative_powell(u_ll, u_rr,
-                                             normal_direction_ll::AbstractVector,
-                                             normal_direction_average::AbstractVector,
+                                             normal_direction::AbstractVector,
                                              equations::IdealGlmMhdEquations2D)
     rho_ll, rho_v1_ll, rho_v2_ll, rho_v3_ll, rho_e_ll, B1_ll, B2_ll, B3_ll, psi_ll = u_ll
     rho_rr, rho_v1_rr, rho_v2_rr, rho_v3_rr, rho_e_rr, B1_rr, B2_rr, B3_rr, psi_rr = u_rr
@@ -265,14 +261,9 @@ end
     v3_ll = rho_v3_ll / rho_ll
     v_dot_B_ll = v1_ll * B1_ll + v2_ll * B2_ll + v3_ll * B3_ll
 
-    # Note that `v_dot_n_ll` uses the `normal_direction_ll` (contravariant vector
-    # at the same node location) while `B_dot_n_rr` uses the averaged normal
-    # direction. The reason for this is that `v_dot_n_ll` depends only on the left
-    # state and multiplies some gradient while `B_dot_n_rr` is used to compute
-    # the divergence of B.
-    v_dot_n_ll = v1_ll * normal_direction_ll[1] + v2_ll * normal_direction_ll[2]
-    B_dot_n_rr = B1_rr * normal_direction_average[1] +
-                 B2_rr * normal_direction_average[2]
+    v_dot_n_ll = v1_ll * normal_direction[1] + v2_ll * normal_direction[2]
+    B_dot_n_rr = B1_rr * normal_direction[1] +
+                 B2_rr * normal_direction[2]
 
     # Powell nonconservative term:   (0, B_1, B_2, B_3, v⋅B, v_1, v_2, v_3, 0)
     # Galilean nonconservative term: (0, 0, 0, 0, ψ v_{1,2}, 0, 0, 0, v_{1,2})
@@ -300,8 +291,9 @@ of the [`IdealGlmMhdEquations2D`](@ref).
 
 This implementation uses a non-conservative term that can be written as the product
 of local and symmetric parts. It is equivalent to the non-conservative flux of Bohm
-et al. (`flux_nonconservative_powell`) for conforming meshes but it yields different
-results on non-conforming meshes(!).
+et al. [`flux_nonconservative_powell`](@ref) for conforming meshes but it yields different
+results on non-conforming meshes(!). On curvilinear meshes this formulation applies the
+local normal direction compared to the averaged one used in [`flux_nonconservative_powell`](@ref).
 
 The two other flux functions with the same name return either the local
 or symmetric portion of the non-conservative flux based on the type of the

src/equations/ideal_glm_mhd_3d.jl

@@ -224,18 +224,15 @@ end
     flux_nonconservative_powell(u_ll, u_rr, orientation::Integer,
                                 equations::IdealGlmMhdEquations3D)
     flux_nonconservative_powell(u_ll, u_rr,
-                                normal_direction_ll     ::AbstractVector,
-                                normal_direction_average::AbstractVector,
+                                normal_direction::AbstractVector,
                                 equations::IdealGlmMhdEquations3D)
 
 Non-symmetric two-point flux discretizing the nonconservative (source) term of
 Powell and the Galilean nonconservative term associated with the GLM multiplier
 of the [`IdealGlmMhdEquations3D`](@ref).
 
-On curvilinear meshes, this nonconservative flux depends on both the
-contravariant vector (normal direction) at the current node and the averaged
-one. This is different from numerical fluxes used to discretize conservative
-terms.
+On curvilinear meshes, the implementation differs from the reference since we use the averaged 
+normal direction for the metrics dealiasing.
 
 ## References
 - Marvin Bohm, Andrew R.Winters, Gregor J. Gassner, Dominik Derigs,
@@ -292,8 +289,7 @@ terms.
 end
 
 @inline function flux_nonconservative_powell(u_ll, u_rr,
-                                             normal_direction_ll::AbstractVector,
-                                             normal_direction_average::AbstractVector,
+                                             normal_direction::AbstractVector,
                                              equations::IdealGlmMhdEquations3D)
     rho_ll, rho_v1_ll, rho_v2_ll, rho_v3_ll, rho_e_ll, B1_ll, B2_ll, B3_ll, psi_ll = u_ll
     rho_rr, rho_v1_rr, rho_v2_rr, rho_v3_rr, rho_e_rr, B1_rr, B2_rr, B3_rr, psi_rr = u_rr
@@ -303,16 +299,11 @@ end
     v3_ll = rho_v3_ll / rho_ll
     v_dot_B_ll = v1_ll * B1_ll + v2_ll * B2_ll + v3_ll * B3_ll
 
-    # Note that `v_dot_n_ll` uses the `normal_direction_ll` (contravariant vector
-    # at the same node location) while `B_dot_n_rr` uses the averaged normal
-    # direction. The reason for this is that `v_dot_n_ll` depends only on the left
-    # state and multiplies some gradient while `B_dot_n_rr` is used to compute
-    # the divergence of B.
-    v_dot_n_ll = v1_ll * normal_direction_ll[1] + v2_ll * normal_direction_ll[2] +
-                 v3_ll * normal_direction_ll[3]
-    B_dot_n_rr = B1_rr * normal_direction_average[1] +
-                 B2_rr * normal_direction_average[2] +
-                 B3_rr * normal_direction_average[3]
+    v_dot_n_ll = v1_ll * normal_direction[1] + v2_ll * normal_direction[2] +
+                 v3_ll * normal_direction[3]
+    B_dot_n_rr = B1_rr * normal_direction[1] +
+                 B2_rr * normal_direction[2] +
+                 B3_rr * normal_direction[3]
 
     # Powell nonconservative term:   (0, B_1, B_2, B_3, v⋅B, v_1, v_2, v_3, 0)
     # Galilean nonconservative term: (0, 0, 0, 0, ψ v_{1,2,3}, 0, 0, 0, v_{1,2,3})

src/equations/shallow_water_2d.jl

@@ -267,8 +267,7 @@ end
     flux_nonconservative_wintermeyer_etal(u_ll, u_rr, orientation::Integer,
                                           equations::ShallowWaterEquations2D)
     flux_nonconservative_wintermeyer_etal(u_ll, u_rr,
-                                          normal_direction_ll     ::AbstractVector,
-                                          normal_direction_average::AbstractVector,
+                                          normal_direction::AbstractVector,
                                           equations::ShallowWaterEquations2D)
 
 Non-symmetric two-point volume flux discretizing the nonconservative (source) term
@@ -303,26 +302,24 @@ Further details are available in the papers:
 end
 
 @inline function flux_nonconservative_wintermeyer_etal(u_ll, u_rr,
-                                                       normal_direction_ll::AbstractVector,
-                                                       normal_direction_average::AbstractVector,
+                                                       normal_direction::AbstractVector,
                                                        equations::ShallowWaterEquations2D)
     # Pull the necessary left and right state information
     h_ll = waterheight(u_ll, equations)
     b_jump = u_rr[4] - u_ll[4]
 
     # Bottom gradient nonconservative term: (0, g h b_x, g h b_y, 0)
     return SVector(0,
-                   normal_direction_average[1] * equations.gravity * h_ll * b_jump,
-                   normal_direction_average[2] * equations.gravity * h_ll * b_jump,
+                   normal_direction[1] * equations.gravity * h_ll * b_jump,
+                   normal_direction[2] * equations.gravity * h_ll * b_jump,
                    0)
 end
 
 """
     flux_nonconservative_fjordholm_etal(u_ll, u_rr, orientation::Integer,
                                         equations::ShallowWaterEquations2D)
     flux_nonconservative_fjordholm_etal(u_ll, u_rr,
-                                        normal_direction_ll     ::AbstractVector,
-                                        normal_direction_average::AbstractVector,
+                                        normal_direction::AbstractVector,
                                         equations::ShallowWaterEquations2D)
 
 Non-symmetric two-point surface flux discretizing the nonconservative (source) term of
@@ -365,8 +362,7 @@ and for curvilinear 2D case in the paper:
 end
 
 @inline function flux_nonconservative_fjordholm_etal(u_ll, u_rr,
-                                                     normal_direction_ll::AbstractVector,
-                                                     normal_direction_average::AbstractVector,
+                                                     normal_direction::AbstractVector,
                                                      equations::ShallowWaterEquations2D)
     # Pull the necessary left and right state information
     h_ll, _, _, b_ll = u_ll
@@ -376,8 +372,8 @@ end
     b_jump = b_rr - b_ll
 
     # Bottom gradient nonconservative term: (0, g h b_x, g h b_y, 0)
-    f2 = normal_direction_average[1] * equations.gravity * h_average * b_jump
-    f3 = normal_direction_average[2] * equations.gravity * h_average * b_jump
+    f2 = normal_direction[1] * equations.gravity * h_average * b_jump
+    f3 = normal_direction[2] * equations.gravity * h_average * b_jump
 
     # First and last equations do not have a nonconservative flux
     f1 = f4 = 0
@@ -424,8 +420,7 @@ end
     flux_nonconservative_audusse_etal(u_ll, u_rr, orientation::Integer,
                                       equations::ShallowWaterEquations2D)
     flux_nonconservative_audusse_etal(u_ll, u_rr,
-                                      normal_direction_ll     ::AbstractVector,
-                                      normal_direction_average::AbstractVector,
+                                      normal_direction::AbstractVector,
                                       equations::ShallowWaterEquations2D)
 
 Non-symmetric two-point surface flux that discretizes the nonconservative (source) term.
@@ -467,8 +462,7 @@ Further details for the hydrostatic reconstruction and its motivation can be fou
 end
 
 @inline function flux_nonconservative_audusse_etal(u_ll, u_rr,
-                                                   normal_direction_ll::AbstractVector,
-                                                   normal_direction_average::AbstractVector,
+                                                   normal_direction::AbstractVector,
                                                    equations::ShallowWaterEquations2D)
     # Pull the water height and bottom topography on the left
     h_ll, _, _, b_ll = u_ll
@@ -479,8 +473,8 @@ end
     # Copy the reconstructed water height for easier to read code
     h_ll_star = u_ll_star[1]
 
-    f2 = normal_direction_average[1] * equations.gravity * (h_ll^2 - h_ll_star^2)
-    f3 = normal_direction_average[2] * equations.gravity * (h_ll^2 - h_ll_star^2)
+    f2 = normal_direction[1] * equations.gravity * (h_ll^2 - h_ll_star^2)
+    f3 = normal_direction[2] * equations.gravity * (h_ll^2 - h_ll_star^2)
 
     # First and last equations do not have a nonconservative flux
     f1 = f4 = 0

src/solvers/dgmulti/dg.jl

@@ -407,11 +407,7 @@ function calc_interface_flux!(cache, surface_integral::SurfaceIntegralWeakForm,
             # Two notes on the use of `flux_nonconservative`:
             # 1. In contrast to other mesh types, only one nonconservative part needs to be
             #    computed since we loop over the elements, not the unique interfaces.
-            # 2. In general, nonconservative fluxes can depend on both the contravariant
-            #    vectors (normal direction) at the current node and the averaged ones. However,
-            #    both are the same at watertight interfaces, so we pass `normal` twice.
-            nonconservative_part = flux_nonconservative(uM, uP, normal, normal,
-                                                        equations)
+            nonconservative_part = flux_nonconservative(uM, uP, normal, equations)
             # The factor 0.5 is necessary for the nonconservative fluxes based on the
             # interpretation of global SBP operators.
             flux_face_values[idM] = (conservative_part + 0.5 * nonconservative_part) *
@@ -556,14 +552,12 @@ function calc_single_boundary_flux!(cache, t, boundary_condition, boundary_key,
                                                         surface_flux, equations)
 
             # Compute pointwise nonconservative numerical flux at the boundary.
-            # In general, nonconservative fluxes can depend on both the contravariant
-            # vectors (normal direction) at the current node and the averaged ones.
-            # However, there is only one `face_normal` at boundaries, which we pass in twice.
-            # Note: This does not set any type of boundary condition for the nonconservative term
-            noncons_flux_at_face_node = nonconservative_flux(u_face_values[i, f],
-                                                             u_face_values[i, f],
-                                                             face_normal, face_normal,
-                                                             equations)
+            noncons_flux_at_face_node = boundary_condition(u_face_values[i, f],
+                                                           face_normal,
+                                                           face_coordinates,
+                                                           t,
+                                                           nonconservative_flux,
+                                                           equations)
 
             flux_face_values[i, f] = (cons_flux_at_face_node +
                                       0.5 * noncons_flux_at_face_node) * Jf[i, f]

src/solvers/dgmulti/flux_differencing.jl

@@ -76,10 +76,7 @@ end
         du_i = du[i]
         for j in col_ids
             u_j = u[j]
-            # The `normal_direction::AbstractVector` has to be passed in twice.
-            # This is because on curved meshes, nonconservative fluxes are
-            # evaluated using both the normal and its average at interfaces.
-            f_ij = volume_flux(u_i, u_j, normal_direction, normal_direction, equations)
+            f_ij = volume_flux(u_i, u_j, normal_direction, equations)
             du_i = du_i + 2 * A[i, j] * f_ij
         end
         du[i] = du_i
@@ -176,11 +173,8 @@ end
         for id in nzrange(A_base, i)
             A_ij = vals[id]
             j = rows[id]
-            # The `normal_direction::AbstractVector` has to be passed in twice.
-            # This is because on curved meshes, nonconservative fluxes are
-            # evaluated using both the normal and its average at interfaces.
             u_j = u[j]
-            f_ij = volume_flux(u_i, u_j, normal_direction, normal_direction, equations)
+            f_ij = volume_flux(u_i, u_j, normal_direction, equations)
             du_i = du_i + 2 * A_ij * f_ij
         end
         du[i] = du_i