update flux_nonconservative_ersing_etal (#53)

Project.toml

@@ -19,5 +19,5 @@ Printf = "1"
 Roots = "2.1.6"
 Static = "0.3, 0.4, 0.5, 0.6, 0.7, 0.8"
 StaticArrays = "1"
-Trixi = "0.7 - 0.8.6"
+Trixi = "0.8.7"
 julia = "1.8"

examples/tree_2d_dgsem/elixir_shallowwater_wall.jl

@@ -30,8 +30,8 @@ boundary_condition = boundary_condition_slip_wall
 ###############################################################################
 # Get the DG approximation space
 
-volume_flux = (flux_wintermeyer_etal, flux_nonconservative_ersing_etal)
-surface_flux = (flux_lax_friedrichs, flux_nonconservative_ersing_etal)
+volume_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal)
+surface_flux = (flux_lax_friedrichs, flux_nonconservative_wintermeyer_etal)
 solver = DGSEM(polydeg = 3, surface_flux = surface_flux,
                volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
 

src/TrixiShallowWater.jl

@@ -25,7 +25,8 @@ export ShallowWaterEquationsWetDry1D, ShallowWaterEquationsWetDry2D,
 
 export hydrostatic_reconstruction_chen_noelle, flux_nonconservative_chen_noelle,
        min_max_speed_chen_noelle, flux_hll_chen_noelle,
-       flux_ersing_etal, flux_es_ersing_etal, hydrostatic_reconstruction_ersing_etal
+       flux_ersing_etal, flux_nonconservative_ersing_etal,
+       flux_es_ersing_etal, hydrostatic_reconstruction_ersing_etal
 
 export ManningFriction, MeyerPeterMueller, GrassModel, ShieldsStressModel,
        dissipation_roe, water_sediment_height, source_term_bottom_friction

src/equations/shallow_water_exner_1d.jl

@@ -347,9 +347,9 @@ inequality.
 This non-conservative flux should be used together with [`flux_ersing_etal`](@ref) to create a
 scheme that is entropy conservative and well-balanced.
 """
-@inline function Trixi.flux_nonconservative_ersing_etal(u_ll, u_rr,
-                                                        orientation::Integer,
-                                                        equations::ShallowWaterExnerEquations1D)
+@inline function flux_nonconservative_ersing_etal(u_ll, u_rr,
+                                                  orientation::Integer,
+                                                  equations::ShallowWaterExnerEquations1D)
     # Pull the necessary left and right state information
     h_ll, _, h_b_ll = u_ll
     h_rr, _, h_b_rr = u_rr

src/equations/shallow_water_multilayer_1d.jl

@@ -310,9 +310,9 @@ In the two-layer setting this combination is equivalent to the fluxes in:
   curvilinear meshes
   [DOI: 10.1007/s10915-024-02451-2](https://doi.org/10.1007/s10915-024-02451-2)
 """
-@inline function Trixi.flux_nonconservative_ersing_etal(u_ll, u_rr,
-                                                        orientation::Integer,
-                                                        equations::ShallowWaterMultiLayerEquations1D)
+@inline function flux_nonconservative_ersing_etal(u_ll, u_rr,
+                                                  orientation::Integer,
+                                                  equations::ShallowWaterMultiLayerEquations1D)
     # Pull the necessary left and right state information
     h_ll = waterheight(u_ll, equations)
     h_rr = waterheight(u_rr, equations)

src/equations/shallow_water_multilayer_2d.jl

@@ -467,9 +467,9 @@ In the two-layer setting this combination is equivalent to the fluxes in:
   curvilinear meshes
   [DOI: 10.1007/s10915-024-02451-2](https://doi.org/10.1007/s10915-024-02451-2)
 """
-@inline function Trixi.flux_nonconservative_ersing_etal(u_ll, u_rr,
-                                                        orientation::Integer,
-                                                        equations::ShallowWaterMultiLayerEquations2D)
+@inline function flux_nonconservative_ersing_etal(u_ll, u_rr,
+                                                  orientation::Integer,
+                                                  equations::ShallowWaterMultiLayerEquations2D)
     # Pull the necessary left and right state information
     h_ll = waterheight(u_ll, equations)
     h_rr = waterheight(u_rr, equations)
@@ -508,10 +508,10 @@ In the two-layer setting this combination is equivalent to the fluxes in:
     return SVector(f)
 end
 
-@inline function Trixi.flux_nonconservative_ersing_etal(u_ll, u_rr,
-                                                        normal_direction_ll::AbstractVector,
-                                                        normal_direction_average::AbstractVector,
-                                                        equations::ShallowWaterMultiLayerEquations2D)
+@inline function flux_nonconservative_ersing_etal(u_ll, u_rr,
+                                                  normal_direction_ll::AbstractVector,
+                                                  normal_direction_average::AbstractVector,
+                                                  equations::ShallowWaterMultiLayerEquations2D)
     # Pull the necessary left and right state information
     h_ll = waterheight(u_ll, equations)
     h_rr = waterheight(u_rr, equations)

src/equations/shallow_water_two_layer_1d.jl

@@ -230,9 +230,9 @@ For further details see:
   curvilinear meshes
   [DOI: 10.1007/s10915-024-02451-2](https://doi.org/10.1007/s10915-024-02451-2)
 """
-@inline function Trixi.flux_nonconservative_ersing_etal(u_ll, u_rr,
-                                                        orientation::Integer,
-                                                        equations::ShallowWaterTwoLayerEquations1D)
+@inline function flux_nonconservative_ersing_etal(u_ll, u_rr,
+                                                  orientation::Integer,
+                                                  equations::ShallowWaterTwoLayerEquations1D)
     # Pull the necessary left and right state information
     h_upper_ll, h_lower_ll = waterheight(u_ll, equations)
     h_upper_rr, h_lower_rr = waterheight(u_rr, equations)

src/equations/shallow_water_two_layer_2d.jl

@@ -343,9 +343,9 @@ For further details see:
   curvilinear meshes
   [DOI: 10.1007/s10915-024-02451-2](https://doi.org/10.1007/s10915-024-02451-2)
 """
-@inline function Trixi.flux_nonconservative_ersing_etal(u_ll, u_rr,
-                                                        orientation::Integer,
-                                                        equations::ShallowWaterTwoLayerEquations2D)
+@inline function flux_nonconservative_ersing_etal(u_ll, u_rr,
+                                                  orientation::Integer,
+                                                  equations::ShallowWaterTwoLayerEquations2D)
     # Pull the necessary left and right state information
     h_upper_ll, h_lower_ll = waterheight(u_ll, equations)
     h_upper_rr, h_lower_rr = waterheight(u_rr, equations)
@@ -381,10 +381,10 @@ For further details see:
     return f
 end
 
-@inline function Trixi.flux_nonconservative_ersing_etal(u_ll, u_rr,
-                                                        normal_direction_ll::AbstractVector,
-                                                        normal_direction_average::AbstractVector,
-                                                        equations::ShallowWaterTwoLayerEquations2D)
+@inline function flux_nonconservative_ersing_etal(u_ll, u_rr,
+                                                  normal_direction_ll::AbstractVector,
+                                                  normal_direction_average::AbstractVector,
+                                                  equations::ShallowWaterTwoLayerEquations2D)
     # Pull the necessary left and right state information
     h_upper_ll, h_lower_ll = waterheight(u_ll, equations)
     h_upper_rr, h_lower_rr = waterheight(u_rr, equations)

src/equations/shallow_water_wet_dry_1d.jl

@@ -214,11 +214,18 @@ end
 Non-symmetric two-point volume flux discretizing the nonconservative (source) term
 that contains the gradient of the bottom topography [`ShallowWaterEquationsWetDry1D`](@ref).
 
-Further details are available in the paper:#include("numerical_fluxes.jl")
+Gives entropy conservation and well-balancedness on both the volume and surface when combined with 
+[`flux_wintermeyer_etal`](@ref).
+
+Further details are available in the papers:
 - Niklas Wintermeyer, Andrew R. Winters, Gregor J. Gassner and David A. Kopriva (2017)
   An entropy stable nodal discontinuous Galerkin method for the two dimensional
   shallow water equations on unstructured curvilinear meshes with discontinuous bathymetry
   [DOI: 10.1016/j.jcp.2017.03.036](https://doi.org/10.1016/j.jcp.2017.03.036)
+- Patrick Ersing, Andrew R. Winters (2023)
+  An entropy stable discontinuous Galerkin method for the two-layer shallow water equations on 
+  curvilinear meshes
+  [DOI: 10.48550/arXiv.2306.12699](https://doi.org/10.48550/arXiv.2306.12699)
 """
 @inline function Trixi.flux_nonconservative_wintermeyer_etal(u_ll, u_rr,
                                                              orientation::Integer,
@@ -234,11 +241,8 @@ end
 Non-symmetric two-point surface flux discretizing the nonconservative (source) term of
 that contains the gradient of the bottom topography [`ShallowWaterEquationsWetDry1D`](@ref).
 
-This contains additional terms compared to [`flux_nonconservative_wintermeyer_etal`](@ref)
-that account for possible discontinuities in the bottom topography function.
-Thus, this flux should be used in general at interfaces. For flux differencing volume terms,
-[`flux_nonconservative_wintermeyer_etal`](@ref) is analytically equivalent but slightly
-cheaper.
+This flux can be used together with [`flux_fjordholm_etal`](@ref) at interfaces to ensure entropy
+conservation and well-balancedness.
 
 Further details for the original finite volume formulation are available in
 - Ulrik S. Fjordholm, Siddhartha Mishr and Eitan Tadmor (2011)
@@ -320,40 +324,12 @@ Further details on the hydrostatic reconstruction and its motivation can be foun
     h_ll_star = u_ll_star[1]
 
     z = zero(eltype(u_ll))
-    # Includes two parts:
-    #   (i)  Diagonal (consistent) term from the volume flux that uses `b_ll` to avoid
-    #        cross-averaging across a discontinuous bottom topography
-    #   (ii) True surface part that uses `h_ll` and `h_ll_star` to handle discontinuous bathymetry
+
     return SVector(z,
-                   equations.gravity * h_ll * b_ll -
-                   equations.gravity * (h_ll_star + h_ll) * (b_ll - b_star),
+                   -equations.gravity * (h_ll_star + h_ll) * (b_ll - b_star),
                    z)
 end
 
-"""
-    flux_nonconservative_ersing_etal(u_ll, u_rr, orientation::Integer,
-                                     equations::ShallowWaterEquationsWetDry1D)
-
-Non-symmetric path-conservative two-point volume flux discretizing the nonconservative (source) term
-that contains the gradient of the bottom topography [`ShallowWaterEquationsWetDry1D`](@ref).
-
-This is a modified version of [`flux_nonconservative_wintermeyer_etal`](@ref) that gives entropy 
-conservation and well-balancedness in both the volume and surface when combined with 
-[`flux_wintermeyer_etal`](@ref).
-
-For further details see:
-- Patrick Ersing, Andrew R. Winters (2023)
-  An entropy stable discontinuous Galerkin method for the two-layer shallow water equations on 
-  curvilinear meshes
-  [DOI: 10.1007/s10915-024-02451-2](https://doi.org/10.1007/s10915-024-02451-2)
-"""
-@inline function Trixi.flux_nonconservative_ersing_etal(u_ll, u_rr,
-                                                        orientation::Integer,
-                                                        equations::ShallowWaterEquationsWetDry1D)
-    return Trixi.flux_nonconservative_ersing_etal(u_ll, u_rr, orientation,
-                                                  equations.basic_swe)
-end
-
 """
     flux_fjordholm_etal(u_ll, u_rr, orientation,
                         equations::ShallowWaterEquationsWetDry1D)
@@ -379,7 +355,8 @@ end
 
 Total energy conservative (mathematical entropy for shallow water equations) split form.
 When the bottom topography is nonzero this scheme will be well-balanced when used as a `volume_flux`.
-The `surface_flux` should still use, e.g., [`flux_fjordholm_etal`](@ref).
+For the `surface_flux` either [`flux_wintermeyer_etal`](@ref) or [`flux_fjordholm_etal`](@ref) can
+be used to ensure well-balancedness and entropy conservation.
 
 Further details are available in Theorem 1 of the paper:
 - Niklas Wintermeyer, Andrew R. Winters, Gregor J. Gassner and David A. Kopriva (2017)

src/equations/shallow_water_wet_dry_2d.jl

@@ -250,16 +250,18 @@ end
 Non-symmetric two-point volume flux discretizing the nonconservative (source) term
 that contains the gradient of the bottom topography [`ShallowWaterEquationsWetDry2D`](@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.
+For the `surface_flux` either [`flux_wintermeyer_etal`](@ref) or [`flux_fjordholm_etal`](@ref) can
+be used to ensure well-balancedness and entropy conservation.
 
-Further details are available in the paper:
+Further details are available in the papers:
 - Niklas Wintermeyer, Andrew R. Winters, Gregor J. Gassner and David A. Kopriva (2017)
   An entropy stable nodal discontinuous Galerkin method for the two dimensional
   shallow water equations on unstructured curvilinear meshes with discontinuous bathymetry
   [DOI: 10.1016/j.jcp.2017.03.036](https://doi.org/10.1016/j.jcp.2017.03.036)
+- Patrick Ersing, Andrew R. Winters (2023)
+  An entropy stable discontinuous Galerkin method for the two-layer shallow water equations on 
+  curvilinear meshes
+  [DOI: 10.48550/arXiv.2306.12699](https://doi.org/10.48550/arXiv.2306.12699)
 """
 @inline function Trixi.flux_nonconservative_wintermeyer_etal(u_ll, u_rr,
                                                              orientation::Integer,
@@ -289,16 +291,8 @@ end
 Non-symmetric two-point surface flux discretizing the nonconservative (source) term of
 that contains the gradient of the bottom topography [`ShallowWaterEquationsWetDry2D`](@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.
-
-This contains additional terms compared to [`flux_nonconservative_wintermeyer_etal`](@ref)
-that account for possible discontinuities in the bottom topography function.
-Thus, this flux should be used in general at interfaces. For flux differencing volume terms,
-[`flux_nonconservative_wintermeyer_etal`](@ref) is analytically equivalent but slightly
-cheaper.
+This flux can be used together with [`flux_fjordholm_etal`](@ref) at interfaces to ensure entropy
+conservation and well-balancedness.
 
 Further details for the original finite volume formulation are available in
 - Ulrik S. Fjordholm, Siddhartha Mishr and Eitan Tadmor (2011)
@@ -487,18 +481,15 @@ Further details on the hydrostatic reconstruction and its motivation can be foun
     h_ll_star = u_ll_star[1]
 
     z = zero(eltype(u_ll))
-    # Includes two parts:
-    #   (i)  Diagonal (consistent) term from the volume flux that uses `b_ll` to avoid
-    #        cross-averaging across a discontinuous bottom topography
-    #   (ii) True surface part that uses `h_ll` and `h_ll_star` to handle discontinuous bathymetry
+
     g = equations.gravity
     if orientation == 1
         f = SVector(z,
-                    g * h_ll * b_ll - g * (h_ll_star + h_ll) * (b_ll - b_star),
+                    -g * (h_ll_star + h_ll) * (b_ll - b_star),
                     z, z)
     else # orientation == 2
         f = SVector(z, z,
-                    g * h_ll * b_ll - g * (h_ll_star + h_ll) * (b_ll - b_star),
+                    -g * (h_ll_star + h_ll) * (b_ll - b_star),
                     z)
     end
 
@@ -524,73 +515,17 @@ end
     # Copy the reconstructed water height for easier to read code
     h_ll_star = u_ll_star[1]
 
-    # Comes in two parts:
-    #   (i)  Diagonal (consistent) term from the volume flux that uses `normal_direction_average`
-    #        but we use `b_ll` to avoid cross-averaging across a discontinuous bottom topography
-
-    f2 = normal_direction_average[1] * equations.gravity * h_ll * b_ll
-    f3 = normal_direction_average[2] * equations.gravity * h_ll * b_ll
-
-    #   (ii) True surface part that uses `normal_direction_ll`, `h_ll` and `h_ll_star`
-    #        to handle discontinuous bathymetry
-
-    f2 -= normal_direction_ll[1] * equations.gravity * (h_ll_star + h_ll) *
-          (b_ll - b_star)
-    f3 -= normal_direction_ll[2] * equations.gravity * (h_ll_star + h_ll) *
-          (b_ll - b_star)
+    f2 = -normal_direction_average[1] * equations.gravity * (h_ll_star + h_ll) *
+         (b_ll - b_star)
+    f3 = -normal_direction_average[2] * equations.gravity * (h_ll_star + h_ll) *
+         (b_ll - b_star)
 
     # First and last equations do not have a nonconservative flux
     f1 = f4 = zero(eltype(u_ll))
 
     return SVector(f1, f2, f3, f4)
 end
 
-"""
-    flux_nonconservative_ersing_etal(u_ll, u_rr, orientation::Integer,
-                                     equations::ShallowWaterEquationsWetDry2D)
-    flux_nonconservative_ersing_etal(u_ll, u_rr,
-                                     normal_direction_ll::AbstractVector,
-                                     normal_direction_average::AbstractVector,
-                                     equations::ShallowWaterEquationsWetDry2D)
-
-!!! warning "Experimental code"
-    This numerical flux is experimental and may change in any future release.
-
-Non-symmetric path-conservative two-point volume flux discretizing the nonconservative (source) term
-that contains the gradient of the bottom topography [`ShallowWaterEquationsWetDry2D`](@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.
-
-This is a modified version of [`flux_nonconservative_wintermeyer_etal`](@ref) that gives entropy 
-conservation and well-balancedness in both the volume and surface when combined with 
-[`flux_wintermeyer_etal`](@ref).
-
-For further details see:
-- Patrick Ersing, Andrew R. Winters (2023)
-  An entropy stable discontinuous Galerkin method for the two-layer shallow water equations on 
-  curvilinear meshes
-  [DOI: 10.1007/s10915-024-02451-2](https://doi.org/10.1007/s10915-024-02451-2)
-"""
-@inline function Trixi.flux_nonconservative_ersing_etal(u_ll, u_rr,
-                                                        orientation::Integer,
-                                                        equations::ShallowWaterEquationsWetDry2D)
-    return Trixi.flux_nonconservative_ersing_etal(u_ll, u_rr, orientation,
-                                                  equations.basic_swe)
-end
-
-@inline function Trixi.flux_nonconservative_ersing_etal(u_ll, u_rr,
-                                                        normal_direction_ll::AbstractVector,
-                                                        normal_direction_average::AbstractVector,
-                                                        equations::ShallowWaterEquationsWetDry2D)
-    Trixi.flux_nonconservative_ersing_etal(u_ll, u_rr,
-                                           normal_direction_ll,
-                                           normal_direction_average,
-                                           equations.basic_swe)
-end
-
 """
     flux_fjordholm_etal(u_ll, u_rr, orientation_or_normal_direction,
                         equations::ShallowWaterEquationsWetDry2D)
@@ -622,7 +557,8 @@ end
 
 Total energy conservative (mathematical entropy for shallow water equations) split form.
 When the bottom topography is nonzero this scheme will be well-balanced when used as a `volume_flux`.
-The `surface_flux` should still use, e.g., [`flux_fjordholm_etal`](@ref).
+For the `surface_flux` either [`flux_wintermeyer_etal`](@ref) or [`flux_fjordholm_etal`](@ref) can
+be used to ensure well-balancedness and entropy conservation.
 
 Further details are available in Theorem 1 of the paper:
 - Niklas Wintermeyer, Andrew R. Winters, Gregor J. Gassner and David A. Kopriva (2017)

test/Project.toml

@@ -12,7 +12,7 @@ OrdinaryDiffEq = "6.49.1"
 Printf = "1"
 Roots = "2.1.6"
 Test = "1"
-Trixi = "0.7 - 0.8.6"
+Trixi = "0.8.7"
 
 [preferences.OrdinaryDiffEq]
 PrecompileAutoSpecialize = false